Hyperbolic compartmental models for epidemic spread on networks with uncertain data: application to the emergence of Covid-19 in Italy

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 $\Omega$ for $t>0$:

	$\displaystyle R_{0}(t)$	$\displaystyle=\frac{\int_{\Omega}f_{I}(S_{T},I_{T})\,dx}{\int_{\Omega}\gamma_{I}(x)I_{T}(x,t)\,dx}\cdot\frac{\int_{\Omega}a(x)\sigma(x)E_{T}(x,t)\,dx}{\int_{\Omega}a(x)E_{T}(x,t)\,dx}$		(2.7)
		$\displaystyle+\frac{\int_{\Omega}f_{A}(S_{T},A_{T})\,dx}{\int_{\Omega}\gamma_{A}(x)A_{T}(x,t)\,dx}\cdot\frac{\int_{\Omega}a(x)(1-\sigma(x))E_{T}(x,t)\,dx}{\int_{\Omega}a(x)E_{T}(x,t)\,dx}\,.$

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$: $\gamma_{I}^{-1}$ and $\gamma_{A}^{-1}$, respectively;

• •

  the fraction of individuals passing from compartment $E$ to $I$ and $A$: $a\sigma/a=\sigma$ and $a(1-\sigma)/a=1-\sigma$, 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}$.

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

$$D_{i}=\lambda_{i}^{2}\tau_{i},\quad i\in\{S,E,I,A,R\}$$

that characterize the diffusive transport mechanism of $S,E,I,A,R$ respectively, and letting $\tau_{i}\to 0$, $i\in\{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

$$\begin{split}J_{S}=-D_{S}\frac{\partial S_{c}}{\partial x},\quad J_{E}=-D_{E}\frac{\partial E_{c}}{\partial x},\quad J_{I}=-D_{I}\frac{\partial I_{c}}{\partial x},\quad J_{A}=-D_{A}\frac{\partial A_{c}}{\partial x},\quad J_{R}=-D_{R}\frac{\partial R_{c}}{\partial x}.\end{split}$$

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

$$\begin{split}\frac{\partial S_{c}}{\partial t}&=\frac{\partial}{\partial x}\left(D_{S}\frac{\partial}{\partial x}S_{c}\right)-f_{I}(S_{c},I_{T})-f_{A}(S_{c},A_{T})\\ \frac{\partial E_{c}}{\partial t}&=\frac{\partial}{\partial x}\left(D_{E}\frac{\partial}{\partial x}E_{c}\right)+f_{I}(S_{c},I_{T})+f_{A}(S_{c},A_{T})-aE_{c}\\ \frac{\partial I_{c}}{\partial t}&=\frac{\partial}{\partial x}\left(D_{I}\frac{\partial}{\partial x}I_{c}\right)+a\sigma E_{c}-\gamma_{I}I_{c}\\ \frac{\partial A_{c}}{\partial t}&=\frac{\partial}{\partial x}\left(D_{A}\frac{\partial}{\partial x}A_{c}\right)+a(1-\sigma)E_{c}-\gamma_{A}A_{c}\\ \frac{\partial R_{c}}{\partial t}&=\frac{\partial}{\partial x}\left(D_{R}\frac{\partial}{\partial x}R_{c}\right)+\gamma_{I}I_{c}+\gamma_{A}A_{c}\,.\end{split}$$

(2.8)

Therefore, the relaxation times can modify the nature of the behavior of the solution [15, 10, 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 $\tau_{i}=\tau_{i}(x)$, $i\in\{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.

A network or a connected graph $\mathcal{G=(N,A)}$ is composed of a finite set of $N$ nodes (or vertices) $\mathcal{N}$ and a finite set of $A$ bidirectional arcs (or edges) $\mathcal{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,\ldots,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\in\mathcal{N}$ having $a_{i}\in\mathcal{A},i\leavevmode\nobreak\ =\leavevmode\nobreak\ 1,\ldots,N_{a,L}$ incoming arcs and $a_{j}\in\mathcal{A},j\leavevmode\nobreak\ =\leavevmode\nobreak\ 1,\ldots,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+\Delta t$, for each of them $\left(1+N_{a,L}\right)$ new states originate at interfaces ($L,i$) and $\left(1+N_{a,0}\right)$ new states originate at interfaces ($0,j$) [47]. To compute them, we need to solve $\left(2+N_{a,L}+N_{a,0}\right)$ 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 $\lambda$ its characteristic velocity, and with $u^{\pm}$ the Riemann Invariants. Introducing transmission coefficients, $\alpha_{i,j}\in[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

$$\begin{split}&u^{*}_{L,i}=u^{+}_{L,i}+\sum_{k=1}^{N_{a,L}}\alpha_{i,k}u_{L,k}^{+}+\alpha_{i,n}u_{n}^{-}\\ &v^{*}_{L,i}=\lambda_{i}\left(u^{+}_{L,i}-\sum_{k=1}^{N_{a,L}}\alpha_{i,k}u_{L,k}^{+}-\alpha_{i,n}u_{n}^{-}\right)\end{split}$$

(2.9)

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

$$\begin{split}&u^{*}_{L,n}=u_{n}^{-}+\sum_{k=1}^{N_{a,L}}\alpha_{n,k}u_{L,k}^{+}+\alpha_{n,n}u_{n}^{-}\\ &v^{*}_{L,n}=-\lambda_{n}\left(u_{n}^{-}-\sum_{k=1}^{N_{a,L}}\alpha_{n,k}u_{L,k}^{+}-\alpha_{n,n}u_{n}^{-}\right).\end{split}$$

(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

$$\begin{split}&u^{*}_{0,j}=u_{0,j}^{-}+\sum_{k=1}^{N_{a,0}}\alpha_{j,k}u_{0,k}^{-}+\alpha_{j,n}u_{n}^{+}\\ &v^{*}_{0,j}=-\lambda_{j}\left(u_{0,j}^{-}-\sum_{k=1}^{N_{a,0}}\alpha_{j,k}u_{0,k}^{-}-\alpha_{j,n}u_{n}^{+}\right).\end{split}$$

(2.11)

and for the node side

$$\begin{split}&u^{*}_{0,n}=u_{n}^{+}+\sum_{k=1}^{N_{a,0}}\alpha_{n,k}u_{0,k}^{-}+\alpha_{n,n}u_{n}^{+}\\ &v^{*}_{0,n}=\lambda_{n}\left(u_{n}^{+}-\sum_{k=1}^{N_{a,0}}\alpha_{n,k}u_{0,k}^{-}-\alpha_{n,n}u_{n}^{+}\right)\end{split}$$

(2.12)

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]

$$v_{L,n}^{*}=\sum_{i=1}^{N_{a,L}}v_{L,i}^{*},\qquad\qquad v_{0,n}^{*}=\sum_{j=1}^{N_{a,0}}v_{0,j}^{*};$$

(2.13)

which is fulfilled imposing at interfaces ($L,i$)

$$\lambda_{i}=\sum_{k=1}^{N_{a,L}}\alpha_{k,i}\lambda_{k}+\alpha_{n,i}\lambda_{n},\qquad\qquad\lambda_{n}=\sum_{k=1}^{N_{a,L}}\alpha_{k,n}\lambda_{k}+\alpha_{n,n}\lambda_{n},$$

(2.14)

and at interfaces ($0,j$)

$$\lambda_{n}=\sum_{k=1}^{N_{a,0}}\alpha_{k,n}\lambda_{k}+\alpha_{n,n}\lambda_{n},\qquad\qquad\lambda_{j}=\sum_{k=1}^{N_{a,0}}\alpha_{k,j}\lambda_{k}+\alpha_{n,j}\lambda_{n}.$$

(2.15)

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,\qquad\qquad u^{*}_{L,n}=u_{n}-\frac{v_{n}}{\lambda_{n}},$$

(2.16)

and at outlet nodes

$$v^{*}_{0,n}=0,\qquad\qquad u^{*}_{0,n}=u_{n}+\frac{v_{n}}{\lambda_{n}}.$$

(2.17)

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\,\mathrm{L}=10^{2}\,\mathrm{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 $\Delta 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: $\Delta x_{n_{1}}=0.025$, $\Delta x_{n_{2}}=0.420$, $\Delta x_{n_{3}}=0.100$, $\Delta x_{n_{4}}=0.085$, $\Delta x_{n_{5}}=0.035$.

Transmission coefficients $\alpha_{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\,\mathrm{T}=1\,\mathrm{day}$, for all the compartments we impose: $\lambda_{a_{1}}=0.4$, $\lambda_{a_{2}}=0.6$, $\lambda_{a_{3}}=0.7$, $\lambda_{a_{4}}=\lambda_{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, $\lambda_{n}=0$, while for the rest of the compartments, namely $S,E,A,R$, $\lambda_{n}=1\,$. In the arcs, the relaxation time is assigned so that the model recovers a hyperbolic regime, hence $\tau_{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 $\tau_{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\sim\mathcal{U}(0,1)$ and the initial conditions for compartment $I$, at each node, are prescribed as

$$I(x,0,z)=i_{0}(1+\mu_{1}z)\,,$$

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 $\mu_{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 [29, 20], we fix $\gamma_{I}=1/14$, $\gamma_{A}=2\gamma_{I}=1/7$, $a=1/3$, considering these clinical parameters deterministic. We also assume the probability rate of developing severe symptoms $\sigma=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 the same compartment, through a simple deterministic SEIAR ODE model set up for the whole Lombardy Region, with the result $\beta_{A}=0.545$. Since after February 23, 2020, Codogno, city in the province of Lodi, is put under strict lockdown as “red zone” [29], a reduced contact rate is considered for the node of Lodi, with $\beta_{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}\approx 10.23\,i_{0}$ and $a_{0}\approx 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}\leavevmode\nobreak\ =\leavevmode\nobreak\ 3.6$, which is in agreement with estimations reported in [29, 20, 55].

On the other hand, the transmission rate of compartment $I$, $\beta_{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:

$$\beta_{I}(0,z)=\beta_{I,0}(1+\mu_{1}\mu_{2}z)\,.$$

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 $\beta_{I,0}=0.03\,\beta_{A}$, as in [29, 20]. Furthermore, $\mu_{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)\,,\qquad A(x,0,z)=9.15\,I(x,0,z)\,,$$

$$S(x,0,z)=N(x)-E(x,0,z)-I(x,0,z)-A(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 $\beta_{A}$ is reduced by 15% (with a consequent update also of $\beta_{I,0}$ and therefore of $\beta_{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, $\beta_{A}$ is ulteriorly reduced by 10% (again with a consequent rearrangement also of $\beta_{I,0}$ and $\beta_{I}$) and coefficients $k$ are increased to $k_{I}=k_{A}=70$.

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 $\mu_{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^{\infty}$ norm is observed in terms of both expected value and variance as the number of collocation points increases.