A viral load-based model for epidemic spread on spatial networks

1.   Introduction

The spreading of infectious diseases is intimately correlated with patterns of human movement, as clearly demonstrated by empirical observations ^[1]. Yet, many compartmental epidemiological models, inspired by the celebrated Kermack-McKendrick theory ^[2], rely on the assumption of homogeneous mixing of the individuals and on the gross number of social contacts among them to provide a big picture of the epidemic trends. Actually, the inclusion of a spatial layer in the mathematical description of epidemic spreading is a quite explored topic in the pertinent literature, cf., e.g., ^[3,4]. Moreover, the recent COVID-19 pandemic has renewed the interest in this issue ^[5,6], given the evident role played by global connections in the worldwide spreading of the epidemic after its localised onset. The spatial layer may be taken into account in essentially two ways. One way consists in introducing a domain, which represents the spatial area of interest, and describing there fluxes of individuals by means of transport-diffusion equations in the spirit of fluid dynamics, cf., e.g., ^[6,7] and [8,Chapter 15]. This approach is generally suited to small spatial scales, because, on one hand, it requires to model local movements of the individuals and, on the other hand, it hardly allows one to account for main directions and preferential paths. An alternative way consists in schematising the spatial structure of the problem by means of a network, whose nodes represent spatial locations and whose links indicate connections among them. Nodes are regarded as spatially homogeneous, hence within them the epidemiological state of the individuals evolves according to contact-and-transmission rules with no spatial structure. However, since individuals are in contact and might also move across the nodes the description is on the whole spatially inhomogeneous and appears to be particularly suited to model main and preferential links. In addition to this, such a spatial representation is quite versatile as it may describe equally well both short-range and long-range connections. Contributions in this direction are e.g., the so-called meta-population models, cf., ^[3,9,10,11]. Here, the movements of the individuals are described as "jumps" from one node of the network to another following the graph edges and the weights prescribed to them, which are interpreted as mobility rates. See also the surveys ^[12,13] and the specific application ^[14] to the case of the COVID-19 epidemics. Instead, we remark that in the aforementioned paper ^[5] the dynamical fluxes of individuals along the links of the network are explicitly modelled in the spirit of the models of vehicular traffic on road networks ^[15,16].

As mentioned at the beginning, another issue in standard compartmental epidemiological models is the general lack of specific disease transmission mechanisms. Contagions are usually modelled only on the basis of the gross numbers of susceptible and infected individuals, assuming that the higher these numbers the higher the probability for a susceptible individual to get in contact with an infected one. This is, for instance, the driving mechanism of the celebrated SIS and SIR models then borrowed by many other more sophisticated compartmental models. Also recent models, partly based on more detailed descriptions of social contacts among the individuals, see e.g., ^[17,18], still rely on a SIR-like structure for the description of the contagion dynamics. Nevertheless, it is known that some specific determinants, such as the viral load carried by the individuals, are at the basis of the infection transmission among the individuals and of its detection through viral tests ^[19]. Measurements of these determinants play also a role in the actuation of confinement measures, such as the quarantine of individuals diagnosed as infected. To the best of our knowledge, epidemiological models describing contagion dynamics by means of a viral load variable are still not common in the literature.

Motivated by the considerations above, in this work we provide as main contribution the design of new mathematical models of contagion dynamics on networks based on viral load transmission. In particular, we propose non-compartmental evolutionary models formalised in terms of the density of individuals in each node of the network and the mean viral load there, which yields a direct measure of the extent of the local contagion. In our models, we follow the idea of the jumps across the nodes, thus we do not describe the fluxes of individuals along the links of the network. To build these models from scratch we adopt a statistical mechanics point of view through Boltzmann-type kinetic equations. This approach allows us to postulate elementary viral load transmission dynamics within the nodes and migration dynamics across the nodes and to upscale them in a formally rigorous way to get the aggregate macroscopic description. We stress that, in our case, a compartmentalisation of the individuals in each node based on their infection condition is not necessary. Indeed, we may rely on the viral load variable to describe contagion dynamics at the individual scale and on the mean viral load to account for the aggregate epidemiological trends. We also remark that kinetic equations on graphs are a still largely under-explored topic in the broad literature on kinetic and statistical models of social interactions, see e.g., ^{[20,21,22,23]}. Therefore, our work provides as a by-product also a contribution in this direction.

In more detail, the paper is organised as follows. In section 2 we introduce the Boltzmann-type kinetic description of social interactions and contagion spread on networks. After obtaining the Boltzmann-type equations accounting for contagion within the nodes and movements across the nodes, we derive exact macroscopic equations for the density of individuals and the mean viral load in each node. We then analyse the large-time trends of these quantities, studying in particular the conditions under which the infection is globally eradicated or blows up depending on the characteristics of the individual contagion dynamics. Next, we specialise the model to the case of commuters, which requires to manage the jumps of the individuals across the nodes in a dedicated manner. In section 3 we extend the previous model by introducing the possibility to put node-wise in isolation (quarantine) individuals diagnosed as infected. For this, we take advantage of the viral load as a descriptor of the microscopic state of the individuals and of the techniques developed in ^[24] for kinetic equations with label switching. In essence, in each node we label the individuals as either non-quarantined or quarantined and we allow for changes of label on the basis of individual viral load levels. Such label switches are described in probability as a result of diagnostic tests, which may or may not detect the infection depending on the viral load carried by the individuals being tested. In section 4 we propose numerical tests of prototypical case studies, which illustrate the type of predictions obtainable from our new models. These tests are also useful to validate the aggregate macroscopic models resulting formally from the kinetic description. Indeed, they show a proper matching between the aggregate trends obtained from the macroscopic models and the original particle dynamics simulated by a Monte Carlo approach. Finally, in section 5 we draw some conclusions and we briefly sketch further possible developments.

2.   Social interactions and contagion spread on networks

2.1.   Boltzmann-type kinetic equations on a graph

Let us consider a large system of interacting individuals migrating on a network (see Figure 1), which we model by a graph with a finite number of vertices and edges. The displacements of the individuals from one vertex to another are described as a Markov-type jump process. Labelling each vertex by an index $i\in \mathcal{I}$, where $\mathcal{I}$ is a finite ordered index set with $\left\lvert{ \mathcal{I}}\right\rvert = n\in\mathbb{N}$, for instance $\mathcal{I} = \{1, \, \dots, \, n\}\subset\mathbb{N}$, we may describe such a stochastic jump process by means of a transition probability

representing the probability to jump from vertex $j$ to vertex $i$. Clearly

hence the square matrix ${\bf{P}}: = [P_{ij}]_{i, j\in \mathcal{I}}\in \mathbb{R}^{n\times n}$, called the transition matrix, is left stochastic as its entries are non-negative and its columns sum to one. Since in our application the nodes of the network represent spatial locations, such as districts, cities, countries depending on the scale of interest, it makes sense to assume that the graph is strongly connected, i.e. it contains a directed path from every vertex to every other vertex. Equivalently, the matrix ${\bf{P}}$ is irreducible.

Besides their current vertex, individuals are further characterised by a microscopic state $v\in \mathbb{R}_+$, which in the application discussed in this paper represents their (dimensionless) viral load understood as a measure of their infectivity. The viral load of the individuals evolves in consequence of social contacts, which, in the spirit of the collisional kinetic theory, we schematise as binary interactions of the form

where $v, \, v_\ast$ represent the pre-interaction viral loads of two interacting individuals and $v', \, v_\ast'$ their post-interaction viral loads. Moreover, $p, \, q\in \mathbb{R}_+$ are either deterministic or stochastic interaction coefficients. The binary rules (2.3) express the contagion dynamics, which we assume to be symmetric. Indeed the second rule is obtained from the first one by switching the roles of $v, \, v_\ast$, which means that there is no preferential order of the individuals in an interaction. We also assume that only individuals in the same vertex of the graph may exchange viral load according to (2.3), whereas individuals in different vertices do not interact. The reason is that contagion dynamics require a certain proximity to be effective. Notice that (2.3) assumes implicitly that such dynamics are the same in each vertex. Different contagion dynamics in different vertices might be taken into account by specifying different interaction coefficients $p_i, \, q_i$, $i\in \mathcal{I}$, from node to node. In particular, we fix

where $\nu_1\in [0, \, 1]$ is the physiological decay rate of the viral load of an infected individual, $\eta$ is a random coefficient modelling stochastic fluctuations of the physiological decay rate and $\nu_2\in [0, \, 1]$ is the transmission rate of the contagion. Denoting by $\langle \cdot \rangle$ the expectation with respect to the law of $\eta$, we assume $\langle \eta \rangle = 0$ and moreover $\eta\geq\nu_1-1$ so as to meet the condition $p\in \mathbb{R}_+$, which, together with $q\in \mathbb{R}_+$, guarantees $v', \, v_\ast'\geq 0$ for all $v, \, v_\ast\geq 0$.

From the particle point of view, we model the microscopic dynamics associated with the processes above as follows. We identify a generic individual by its location $X_t\in \mathcal{I}$ in the network and its viral load $V_t\in \mathbb{R}_+$ at time $t$. Then we describe the evolution of $X_t$, $V_t$ during a time step $\Delta{t} > 0$ by means of the following update rules:

where:

i) $\Theta, \, \Xi\in\{0, \, 1\}$ are independent Bernoulli random variables discriminating whether a vertex jump or a change in the viral load take place ($\Theta, \, \Xi = 1$) or not ($\Theta, \, \Xi = 0$) in the time step $\Delta{t}$. In particular, we let

$\chi, \, \mu > 0$ being the mobility and social contact rates, and we assume $\Delta{t}\leq\min\{\frac{1}{\chi}, \, \frac{1}{\mu}\}$ for consistency;

ii) $J_t\in \mathcal{I}$ is a random variable returning the new vertex after a vertex jump, with

iii) $V'_t\in \mathbb{R}_+$ is the new viral load after a social contact. In view of (2.3) we have $V'_t = pV_t+qV^\ast_t$, where $V^\ast_t$ denotes the viral load at time $t$ of the other individual participating in the interaction;

iv) $\delta_{X_t, X^\ast_t}$ is the Kronecker delta:

which we use to express the fact that only individuals in the same vertex may interact and exchange viral load. Here, $X^\ast_t\in \mathcal{I}$ is the vertex occupied at time $t$ by the other individual involved in the interaction.

To pass from the particle model (2.5) to an aggregate description, which is more suited to the investigation of the emerging collective trends, we introduce the distribution function of the viral load $v\in \mathbb{R}_+$ at time $t > 0$ in vertex $i\in \mathcal{I}$, that we denote by $f_i = f_i(v, t): \mathbb{R}_+\times \mathbb{R}_+\to \mathbb{R}_+$. Taking advantage of the procedure detailed in ^[24], see also ^[25] for further reference, in the continuous-time limit $\Delta{t}\to 0^+$ we find that $f_i$ satisfies formally the following Boltzmann-type kinetic equation in weak form:

where $\varphi: \mathbb{R}_+\to \mathbb{R}$ is an arbitrary observable quantity (test function) and $Q(f_i)$ is the so-called collisional operator (in the jargon of classical kinetic theory), which here describes the effect of binary interactions among pairs of individuals belonging to the same vertex of the graph. More specifically:

where $v'$ is given by (2.3)–(2.4) as a function of the integration variables $v$, $v_\ast$.

2.2.   Aggregate description of mobility and contagion

We observe that (2.7) is a non-conservative kinetic equation, indeed the density of individuals in a vertex of the graph, i.e.

is in general not conserved in time because of the jumps across the vertices. Setting $\varphi(v) = 1$ in (2.7) we obtain that $\rho_i$ satisfies the equation

which in vector notation reads

with $\boldsymbol{\rho}: = (\rho_i)_{i\in \mathcal{I}}$. From (2.10) we can investigate the stationary mass distribution $\boldsymbol{\rho}^\infty\in \mathbb{R}^n_+$ emerging for large times, namely the vector satisfying the equation

In practice, $\boldsymbol{\rho}^\infty$ should be an eigenvector of ${\bf{P}}$ corresponding to the eigenvalue $1$. Notice that ${\bf{P}}$ admits indeed this eigenvalue, which coincides also with its spectral radius, because it is a stochastic matrix. Moreover, since ${\bf{P}}$ is irreducible from Perron-Frobenius theory it follows that the eigenvalue $1$ is simple and that there exists a corresponding eigenvector with strictly positive components. See e.g., ^[26] for details. Such an eigenvector is our candidate for $\boldsymbol{\rho}^\infty$, but we need to manage its non-uniqueness (scalar multiples are also eigenvectors). For this, we notice that in our case the $\ell^1$-norm of $\boldsymbol{\rho}$, i.e. the quantity $\| \boldsymbol{\rho}\|_1: = \sum_{i\in \mathcal{I}}\rho_i$, is conserved in time, indeed from (2.9) and (2.2) we deduce

Thus $\| \boldsymbol{\rho}^\infty\|_1$ is fixed by the initial condition, which removes the remaining degree of freedom in the identification of a unique physically admissible $\boldsymbol{\rho}^\infty\in \mathbb{R}^n_+$ solving (2.11).

Still from Perron-Frobenius theory we have that $1$ is also the spectral limit of ${\bf{P}}$, namely the maximum real part of its eigenvalues. Therefore the real parts of the eigenvalues of the matrix ${\bf{P}}- {\bf{I}}$ are non-positive and, in particular, the unique eigenvalue with null real part is simple. This implies that $\boldsymbol{\rho}^\infty$ is a stable equilibrium of (2.10). Actually, it is also attractive because the eigenvalue with null real part is simply associated with the conservation of the $\ell^1$-norm discussed above.

Summarising, we have proved:

Proposition 2.1. There exists a unique physically admissible solution $\boldsymbol{\rho}^\infty\in \mathbb{R}^n_+$ to $(2.11)$, which is a stable and attractive asymptotic density distribution for $(2.10)$.

Conversely, setting $\varphi(v) = v$ in (2.7) we may investigate the evolution of the mean viral load

which, recalling also (2.3)–(2.4), turns out to satisfy the equation

From here, invoking (2.2) we discover

which we may use as a basis to study the large time trend of the $m_i$'s. In particular, we are interested in the stability of the asymptotic state ${\bf{m}}^\infty = {\bf{0}}$, which represents the eradication of the infection in all nodes of the network. Linearising the previous equation around the equilibrium configuration $(\boldsymbol{\rho}, \, {\bf{m}}) = (\boldsymbol{\rho}^\infty, \, {\bf{0}})$ we obtain

Next, recalling that $\boldsymbol{\rho}^\infty$ has strictly positive components we set $c: = \min_{i\in \mathcal{I}}\rho_i^\infty > 0$ and we use

to conclude that:

i) if $\nu_1 > \nu_2$ then

which for $t\to+\infty$ implies $\sum_{i\in \mathcal{I}}\rho_i^\infty m_i\to 0$ and therefore $m_i\to 0$ for all $i\in \mathcal{I}$. In this case the configuration ${\bf{m}}^\infty = {\bf{0}}$ is locally asymptotically stable, hence the infection may be globally eradicated in the long run;

ii) if $\nu_1 < \nu_2$ then

which for $t\to+\infty$ yields $\sum_{i\in \mathcal{I}}\rho_i^\infty m_i\to+\infty$ and therefore $m_i\to+\infty$ for some $i\in \mathcal{I}$. In this case the configuration ${\bf{m}}^\infty = {\bf{0}}$ is unstable: the infection is not under control and tends to grow in time;

iii) if $\nu_1 = \nu_2$ then $\sum_{i\in \mathcal{I}}\rho_im_i$ is constant in time. This indicates that the configuration ${\bf{m}}^\infty = {\bf{0}}$ is stable but not attractive: the infection is under control but cannot be eradicated and becomes endemic.

2.2.1.   Basic reproduction number $\mathcal{R}_0$

A standard parameter normally invoked in epidemiological models to assess the conditions for an epidemic outbreak is the basic reproduction number $\mathcal{R}_0\geq 0$. It represents the mean number of secondary infections caused by a single infected individual in a population of fully susceptible individuals. From the mathematical point of view, this parameter is associated with the (linear) stability (if $\mathcal{R}_0 < 1$) or instability (if $\mathcal{R}_0 > 1$) of the disease-free equilibrium, i.e. the one in which no infected individuals are present in the population.

The concept of basic reproduction number is typical of compartmental population dynamics models, while for our non-compartmental viral load-based model its definition might be less straightforward. Here we propose a conceptual approach leading to a possible definition of $\mathcal{R}_0$ in our context, which relies on the idea of stability/instability of the asymptotic state ${\bf{m}}^\infty = {\bf{0}}$.

From (2.12) together with (2.9) we deduce the following equations for the mean viral loads $m_i$:

which we now linearise around the equilibrium $(\boldsymbol{\rho}, \, {\bf{m}}) = (\boldsymbol{\rho}^\infty, \, {\bf{0}})$. Writing

where $\epsilon > 0$ is a small parameter, and plugging into (2.14) we obtain, at the leading order in $\epsilon$, the following equations for the perturbations $\tilde{m}_i$:

Introducing the diagonal matrix ${\bf{R}}: = \operatorname{diag}(\rho^\infty_1, \, \dots, \, \rho^\infty_n)$, this linear system may be rewritten in compact form as

where $\tilde{ {\bf{m}}}: = (\tilde{m}_1, \, \dots, \, \tilde{m}_n)$, whence we deduce that the stability of the asymptotic state ${\bf{m}}^\infty = {\bf{0}}$ depends on the spectral properties of the matrix

Notice that ${\bf{A}}$ has the form ${\bf{A}} = {\bf{B}}- {\bf{D}}$ with

where ${\bf{B}}$, ${\bf{D}}$ are both non-negative and moreover ${\bf{D}}$ is diagonal and invertible (at least for either $\chi > 0$ or $\mu\nu_1 > 0$). From Perron-Frobenius theory it is known that $\tilde{ {\bf{m}}} = {\bf{0}}$ is a stable equilibrium of (2.15), and hence ${\bf{m}}^\infty = {\bf{0}}$ is a (linearly) stable equilibrium of (2.14), if and only if the spectral radius of the matrix ${\bf{B}} {\bf{D}}^{-1}$ is smaller than $1$. Conversely, if such a spectral radius is larger than $1$ then $\tilde{ {\bf{m}}} = {\bf{0}}$ and ${\bf{m}}^\infty = {\bf{0}}$ are unstable. We conclude therefore that the spectral radius of ${\bf{B}} {\bf{D}}^{-1}$ may appropriately play the role of $\mathcal{R}_0$ in our context.

Let $\beta_{ij}\in \mathbb{R}$ be the $ij$-entry of the matrix ${\bf{B}} {\bf{D}}^{-1}$, then:

According to Perron-Frobenius theorem,

thus we may estimate $\mathcal{R}_0$ by essentially computing $\sum_{j = 1}^{n}\beta_{ij}$. To this purpose, we introduce the following quantities:

and we estimate:

where we have used that $\sum_{j = 1}^{n}P_{ij}\rho^\infty_j = \rho^\infty_i$, cf. (2.11). Similarly, we also estimate:

so that finally from (2.16) we conclude

Estimates (2.17) may be equivalently rewritten as

whence we observe that they imply either $\mathcal{R}_0 < 1$ or $\mathcal{R}_0 > 1$ depending on whether

respectively. Since $\underline{\rho}^\infty/ \overline{\rho}^\infty\leq 1$, this result is less sharp than the one we have obtained in section 2.2, meaning that the estimates (2.17) of $\mathcal{R}_0$ are quantitatively suboptimal. Nevertheless, they are qualitatively informative as far as the dependence of $\mathcal{R}_0$ on the microscopic parameters of the model is concerned. Notice moreover that if $\underline{\rho}^\infty = \overline{\rho}^\infty$, i.e. if the asymptotic density distribution is the uniform one across the network, then from (2.17) we recover $\nu_2 < \nu_1\Rightarrow \mathcal{R}_0 < 1$ and similarly $\nu_2 > \nu_1\Rightarrow \mathcal{R}_0 > 1$. If, in addition, the network reduces to independent isolated nodes, i.e. if $\chi = 0$ so that individuals do not migrate across the nodes, then (2.17) returns $\mathcal{R}_0 = \nu_2/\nu_1$.

2.2.2.   Node-dependent decay and transmission rates and hydrodynamic limit on the network

If the decay and transmission rates $\nu_1$, $\nu_2$ are not constant but vary from vertex to vertex of the graph, i.e. $\nu_1 = \nu_{1, i}$ and $\nu_2 = \nu_{2, i}$, then (2.13) modifies as

Apart from the relatively trivial cases in which either $\nu_{1, i} > \nu_{2, i}$, $\nu_{1, i} < \nu_{2, i}$ or $\nu_{1, i} = \nu_{2, i}$ for all $i\in \mathcal{I}$, which actually reproduce the already considered scenarios with constant rates, this equation is not particularly informative about the large time trend of the infection. We may get more detailed information by looking instead at the hydrodynamic regime of the kinetic equation (2.7), namely the one in which local interactions within the vertices of the graph are much more frequent than jumps from node to node. This amounts to scaling

where $0 < \epsilon\ll 1$ is a small parameter playing conceptually the role of the Knudsen number of classical kinetic theory. In the hydrodynamic limit $\epsilon\to 0^+$ a splitting of (2.7) is possible in:

i) intra-vertex interactions without jumps on the quick time scale $\tau: = t/\epsilon$, ruled by

which lead $f_i$ to converge for large $\tau$ to a local equilibrium distribution $M_i$ (the local Maxwellian) parametrised by the macroscopic quantities conserved by the interactions (the collisional invariants) in the $i$th vertex;

ii) jumps across the vertices without interactions on the slow time scale $t$, ruled by

which determine the macroscopic evolution of the quantities conserved by the local intra-vertex interactions.

Taking $\varphi(v) = 1$ in (2.18) we discover that in all vertices the density $\rho_i$ is constant on the time scale $\tau$. Taking instead $\varphi(v) = v$ we get

Let us assume first that $\nu_{1, i} > \nu_{2, i}$ for all $i\in \mathcal{I}\setminus\{i^\ast\}$ while $\nu_{1, i^\ast} = \nu_{2, i^\ast}$. Then (2.20) implies $m_i\to 0$ exponentially fast for $\tau\to+\infty$ if $i\neq i^\ast$ while $m_{i^\ast}$ is constant on the $\tau$-scale. Thus for $i\neq i^\ast$ the local Maxwellian is $M_{\rho_i}(v) = \rho_i\delta(v)$, $\delta$ denoting the Dirac delta centred in $v = 0$, whereas for $i = i^\ast$ the local Maxwellian is a distribution $M_{\rho_{i^\ast}, m_{i^\ast}}(v)$ parametrised by both $\rho_{i^\ast}$ and $m_{i^\ast}$. Plugging these local Maxwellians into (2.19) and choosing $\varphi(v) = 1$ yields the same time evolution of the densities $\rho_i$, $i\in \mathcal{I}$, as (2.9). Then, in particular, $\rho_i\to\rho^\infty_i > 0$ for all $i\in \mathcal{I}$ when $t\to +\infty$ as asserted by Proposition 2.1. Choosing instead $\varphi(v) = v$ and $i = i^\ast$ we may investigate the evolution of the mean viral load in the $i^\ast$th vertex, the only one which remains constant during the intra-vertex interactions:

Since the graph is strongly connected we have in particular $P_{i^\ast i^\ast} < 1$, whence $\rho_{i^\ast}m_{i^\ast}\to 0$ for $t\to+\infty$. Considering that $\rho_{i^\ast}$ converges to a non-zero asymptotic value, we finally infer $m_{i^\ast}\to 0$. Thus, in the long run, thanks to the mobility of the individuals across the nodes of the network the infection is eradicated also in the node where local contagion dynamics balance with the physiological recovery. Notice however that such an eradication is slower than in the other nodes, indeed it happens on the $t$-time scale rather than on the quicker $\tau$-scale.

Next, let us consider the opposite situation, namely $\nu_{1, i} = \nu_{2, i}$ for all $i\in \mathcal{I}\setminus\{i^\ast\}$ and $\nu_{1, i^\ast} > \nu_{2, i^\ast}$. In this case, in all vertices $i\neq i^\ast$ the local contagion dynamics yield local Maxwellians $M_{\rho_i, m_i}(v)$ parametrised by both $\rho_i$ and $m_i$ while in the $i^\ast$th vertex it results $m_{i^\ast}\to 0$ as $\tau\to+\infty$, thus the local Maxwellian is $M_{\rho_{i^\ast}}(v) = \rho_{i^\ast}\delta(v)$. Consequently from (2.19) with $\varphi(v) = v$ we deduce

whence, summing both sides over $i\in \mathcal{I}\setminus\{i^\ast\}$,

where we have used the fact that $\sum_{i\in \mathcal{I}\setminus\{i^\ast\}}P_{ij} = 1-P_{i^\ast j}$. Let us introduce the set of indices of the vertices different from vertex $i^\ast$ and directly linked to the latter:

Notice that $\mathcal{I}^\ast$ is non-empty due to the strong connectivity of the graph (if $P_{i^\ast j} = 0$ for all $j\neq i^\ast$ then vertex $i^\ast$ could not be reached from any other vertex). Let $a: = \min_{j\in \mathcal{I}^\ast}P_{i^\ast j} > 0$, then from (2.21) we deduce

whence, integrating both sides in time and considering that $\sum_{i\in \mathcal{I}\setminus\{i^\ast\}}\rho_im_i\geq \sum_{i\in \mathcal{I}^\ast}\rho_im_i$ because $\mathcal{I}^\ast\subseteq \mathcal{I}\setminus\{i^\ast\}$,

By Grönwall's inequality we deduce then

which says that $\sum_{i\in \mathcal{I}^\ast}\rho_im_i\to 0$, hence that $m_i\to 0$ for all $i\in \mathcal{I}^\ast$ as $t\to+\infty$. We conclude that, owing to the mobility of the individuals, the infection is certainly eradicated at least in the vertices directly linked to the one where contagion dynamics are weaker than the physiological recovery of the individuals. In particular, if $\mathcal{I}^\ast = \mathcal{I}\setminus\{i^\ast\}$, i.e. if all vertices of the graph are directly linked to vertex $i^\ast$, then the infection is eradicated in the whole network. It is worth stressing that instead in classical homogeneous kinetic models featuring binary interactions of the form (2.3) with coefficients given by (2.4) there is no way to obtain a convergence to zero of the mean if $\nu_1 = \nu_2$.

Finally, if $\nu_{1, i^\ast} < \nu_{2, i^\ast}$ for some $i^\ast\in \mathcal{I}$ then from (2.20) we easily compute that the intra-vertex interactions produce locally a blow-up of the mean viral load: $m_{i^\ast}\to+\infty$ as $\tau\to+\infty$. The mean viral load will consequently blow in every other vertex $i\in \mathcal{I}$ where $\nu_{1, i}\leq\nu_{2, i}$ while the infection will be eradicated only in the vertices where $\nu_{1, i} > \nu_{2, i}$.

2.3.   The problem of commuters

A particularly interesting variation on the model introduced in the previous sections is the case of commuters, namely individuals who do not travel generically between any two vertices of the graph but mainly between two specific vertices representing e.g., the places where they live and work. For the sake of simplicity, we will assume that they travel only between those two vertices.

A possible particle description of the problem is as follows. Let $(X_{o, t}, \, X_{d, t})\in \mathcal{I}\times \mathcal{I}$ be the origin-destination pair of vertices of a commuter who at time $t$ is in $X_{o, t}$ and is heading for $X_{d, t}$. Let moreover $V_t\in \mathbb{R}_+$ be their viral load. We model the evolution of $(X_{o, t}, \, X_{d, t})$ and $V_t$ during a time step $\Delta{t} > 0$ as

where $\Theta$, $\Xi$ are like in (2.6), $V'_t = pV_t+qV^\ast_t$ from (2.3) and $V^\ast_t$, $X^\ast_{o, t}$ are the viral load and the origin vertex of the other individual participating in the interaction. The vector $(X'_{o, t}, \, X'_{d, t})\in \mathcal{I}\times \mathcal{I}$ expresses the new origin-destination pair of vertices of the commuter after a possible journey from $X_{o, t}$ to $X_{d, t}$. Notice that, by definition of commuter, we may only have either $X'_{o, t} = X_{d, t}$ and $X'_{d, t} = X_{o, t}$ if the journey takes place or $X'_{o, t} = X_{o, t}$ and $X'_{d, t} = X_{d, t}$ if it does not. Therefore, if we denote by $p_{X_{d, t}, X_{o, t}}\in [0, \, 1]$ the probability for such a commuter to actually travel from $X_{o, t}$ to $X_{d, t}$ then the law of the random vector $(X'_{o, t}, \, X'_{d, t})$ is

We are now in a position to deduce from (2.22) an aggregate description in terms of the distribution function of the viral load $v\in \mathbb{R}_+$ of commuters between vertices $i, \, j\in \mathcal{I}$ (in the following we will refer to them as $ij$-commuters for brevity) who at time $t > 0$ are in vertex $i$ and are heading for vertex $j$. Let us denote by $f_{ji} = f_{ji}(v, t): \mathbb{R}_+\times \mathbb{R}_+\to \mathbb{R}_+$ such a distribution function. Using again the technique illustrated in ^[24] with minimal adaptations, in the continuous-time limit $\Delta{t}\to 0^+$ we formally obtain from (2.22) that $f_{ji}$ satisfies the following Boltzmann-type kinetic equation in weak form:

where $v'$ is given by (2.3)–(2.4) and

is the cumulative distribution function of the viral load of the individuals who at time $t$ are in vertex $i$ (independently of the other vertex to which they commute).

Equation (2.23) features two important differences with respect to (2.7), that we now elucidate.

i) The commuting term

allows only for exchanges between $ij$-commuters in $j$ heading for $i$ ($p_{ij}f_{ij}$, gain) and $ij$-commuters in $i$ heading for $j$ ($p_{ji}f_{ji}$, loss). Notice indeed that, by definition of $ij$-commuters, individuals arriving in vertex $i$ from vertices different from $j$ do not contribute to the time variation of $f_{ji}$.

ii) The contagion term

takes into account binary interactions between $ij$-commuters in $i$ and all other individuals in vertex $i$, independently of the vertex to which they commute, because they all contribute equally to the spread of the contagion in vertex $i$.

Remark 2.2. We stress that, in general, $p_{ij}\neq P_{ij}$, cf., (2.1). Indeed, in the model introduced in section 2.1, $P_{ij}$ is the probability that any individual in vertex $j$ jumps to vertex $i$ while here $p_{ij}$ is the probability that an $ij$-commuter currently in vertex $j$ jumps to vertex $i$ (all individuals in vertex $j$ might not be $ij$-commuters).

If we define by

the density of $ij$-commuters in vertex $i$ at time $t$ then from (2.23) with $\varphi(v) = 1$ we obtain the macroscopic equation of the mobility on the graph:

Notice that the time evolution of $\rho_{ji}$ depends only on $\rho_{ji}$ itself and on $\rho_{ij}$, which is the density of $ij$-commuters travelling in the opposite direction. This is consistent with the definition of commuters. Moreover, from the analogous equation for $\rho_{ij}$ we easily see that the total density of $ij$-commuters, namely $\rho_{ji}+\rho_{ij}$, is constant in time. Let us denote by^* $\bar{\rho}_{\{i, j\}}: = \rho_{ji}+\rho_{ij}$ the total mass of $ij$-commuters. Then $\rho_{ij} = \bar{\rho}_{\{i, j\}}-\rho_{ji}$, whence

^*We use the set-style subscript $\{i, j\}$ to mean that, unlike $f_{ji}$ and $\rho_{ji}$, in the notation $\bar{\rho}_{\{i, j\}}$ there is no origin-destination ordering of the indices $i$, $j$ because $\bar{\rho}_{\{i, j\}}$ is the total mass of $ij$-commuters independently of the direction of their journey.

which implies that $\rho_{ji}$, $\rho_{ij}$ reach exponentially fast the asymptotic values

From this, summing over all destination vertices, we may also compute the asymptotic value of

i.e. the total density of individuals in vertex $i$ at time $t$. We have:

Remark 2.3. From (2.24), the evolution equation for the total density $\rho_i$ turns out to be

which is here the counterpart of (2.9). If we assume that all commuters always travel from their origin to their destination, i.e. that no commuter might remain in the origin vertex without travelling, then we have $p_{ji} = 1$ for all $i, \, j\in \mathcal{I}$ and this equation becomes

which reminds more closely of (2.9) but with $P_{ij}\rho_j$ replaced by $\rho_{ij}$. Notice that, unlike (2.9), this is not a self-consistent equation in terms of the $\rho_i$'s alone because commuting requires to keep track of origin and destination vertices.

Likewise, if we define by

the mean viral load of $ij$-commuters in vertex $i$ at time $t$ then from (2.23) with $\varphi(v) = v$ and taking also (2.3)–(2.4) and (2.24) into account we obtain

where

is the total mean viral load in vertex $i$. This term couples macroscopically the contagion dynamics of $ij$-commuters to those of all commuters in $i$ travelling to other vertices different from $j$. This reflects the fact that contagion dynamics are not confined to the transport of viral load along the commuting routes (cf. the first term on the right-hand side of (2.26)). Contagion may spread along all routes because of mixed social contacts in each vertex.

3.   Social interactions with quarantine in the nodes

We consider now the case in which individuals may be quarantined in each node of the network if they are diagnosed as infected. Quarantined individuals do not have social contacts with other individuals, therefore they only diminish their viral load and are eventually readmitted in the society when they are diagnosed as recovered by a subsequent viral test.

To approach this problem we take advantage of the framework introduced in ^[24], where kinetic equations with label switching are discussed. In essence, in each vertex $i\in \mathcal{I}$ of the graph we label the individuals by $h = 1$ if they are not quarantined and by $h = 2$ if instead they are. Next, we introduce the label switch probability

namely the probability that an individual in vertex $i$ with viral load $v$ changes their label from $k$ to $h$. In particular, the label switching $1\to 2$ means that the individual is quarantined; conversely, the label switching $2\to 1$ means that the individual is released from quarantine. When $h = k$ the individual does not change their quarantine state. These probabilities fulfil

From the particle point of view, we now characterise the state of a representative individual at time $t$ by the vertex $X_t\in \mathcal{I}$ to which they belong, the label $H_t\in\{1, \, 2\}$ identifying their quarantine state and their viral load $V_t\in \mathbb{R}_+$. Next, we model the state update in a time step $\Delta{t} > 0$ as

where $\Lambda\in\{0, \, 1\}$ is a Bernoulli random variable such that

which discriminates whether a label switching takes place ($\Lambda = 1$) or not ($\Lambda = 0$) in the time step $\Delta{t}$ and $\lambda > 0$ is the rate of viral testing. Moreover, $I_t\in\{1, \, 2\}$ is the new label assigned to the individual after a viral test, with in particular

The random variables $\Theta$, $\Xi$ are like in section 2.1 but in this case we fix $\Delta{t}\leq\min\{\frac{1}{\chi}, \, \frac{1}{\lambda}, \, \frac{1}{\mu}\}$ for consistency. A remarkable difference with respect to the model in section 2.1 is instead the structure of the random variables $J_t$, $V'_t$ for now only a non-quarantined individual ($h = 1$) can jump to another vertex and can experience social contacts producing viral load changes according to formula (2.3). To take these facts into account we modify the law of $J_t$ as follows:

where the $P^{1}_{ij}$'s coincide with the $P_{ij}$'s of section 2.1 while

because a quarantined individual cannot change vertex. Parallelly, for quarantined ($h = 2$) individuals we convert the binary interaction rule (2.3)–(2.4) into a rule expressing a progressive reduction of viral load due to the lack of social contacts:

which is of the form (2.3) with $p = 1-\gamma+\xi$ and $q = 0$. Here, $\gamma\in [0, \, 1]$ is the decay rate of the viral load of a quarantined individual, which in general is expected to satisfy $\gamma\geq\nu_1$ because quarantined individuals may be exposed to specific medical treatments in addition to the physiological recovery. Moreover, $\xi\in \mathbb{R}$ is a centred random coefficient modelling stochastic fluctuations of the decay rate. With the restriction $\xi\geq\gamma-1$ we guarantee $v'\geq 0$ for all $v\geq 0$. On the whole, we set

where $H^\ast_t\in\{1, \, 2\}$ denotes the label identifying the quarantine state of the other individual participating in the interaction.

Let now $f^h_i = f^h_i(v, t): \mathbb{R}_+\times \mathbb{R}_+\to \mathbb{R}_+$ be the distribution function of the viral load $v\in \mathbb{R}_+$ at time $t > 0$ in the vertex $i\in \mathcal{I}$ for an individual with label $h\in\{1, \, 2\}$. Combining (2.7) and the results in ^[24] we deduce from (3.2) that in the continuous-time limit $\Delta{t}\to 0^+$ the $f^h_i$'s satisfy formally the following kinetic equation in weak form:

where the collisional operator $Q^h$ is labelled by $h$ to reflect the different rules of change of viral load followed by non-quarantined and quarantined individuals. Specifically, the collisional operator $Q^1$ coincides with the bilinear one given in (2.8) whereas $Q^2$ takes the form of a linear collision operator, i.e. it depends on the term $f^h_i(v, t)$ alone rather than on the product $f^h_i(v, t)f^h_i(v_\ast, t)$. Explicitly, the equations for the two compartments $h = 1, \, 2$ read then:

● non-quarantined individuals ($h = 1$):

with $v'$ given by (2.3)–(2.4);

● quarantined individuals ($h = 2$):

with $v'$ given by (3.3).

3.1.   Hydrodynamic model with quarantine on the network

We observe that (3.4) is a fully non-conservative kinetic equation, indeed it does not even conserve the mass of the individuals belonging to each compartment. Let

be the density of individuals in vertex $i$ and compartment $h$ at time $t$. Setting $\varphi(v) = 1$ in (3.4) and observing that $\int_{ \mathbb{R}_+}Q^h(f^h_i)(v, t)\, dv = 0$ for $h = 1, \, 2$ because both microscopic rules (2.3)–(2.4) and (3.3) conserve the number of individuals, we obtain that $\rho^h_i$ satisfies the equation

which depends explicitly on the intra-vertex label switching and the jumps from vertex to vertex.

Remark 3.1. Summing (3.7) over $h$ and defining the global density of individuals in vertex $i$ at time $t$ as

we find

where we have recalled (3.1). If the transition probabilities from vertex to vertex do not depend on $h$ then this equation coincides with (2.9). Otherwise it describes different aggregate dynamics, that cannot be expressed in terms of the $\rho_i$'s alone. This is the case of the model with quarantine, as we have in general $P^1_{ij}\neq P^2_{ij}$ due to different mobility features of the individuals in the two compartments $h = 1, \, 2$.

Remark 3.1 implies that (3.7) is the reference equation for the evolution of the density of individuals on the network in the model with quarantine. Notice that it requires, in general, the knowledge of the kinetic distribution functions $f^k_i$ to compute the label switching at the macroscopic scale. Nevertheless, in the particular case that the label switch probabilities $T^{hk}_i$ are constant with respect to $v$, namely if we assume that the probability to diagnose an individual as either infected or recovered may be considered approximately independent of their viral load, (3.7) becomes instead a self-consistent equation for the densities $\rho^h_i$:

which for $h = 1, \, 2$ and recalling (3.1) produces the system

System (3.8) provides self-consistent macroscopic information on the mass distribution of the individuals on the network but not specifically on the aggregate trends of the infection itself. To get a more complete picture of the epidemic spread, we may use (3.5)–(3.6) to deduce a coupled system of hydrodynamic equations for the mean viral loads

Letting $\varphi(v) = v$ in (3.5)–(3.6) and assuming again that the $T^{hk}_i$'s are constant with respect to $v$ we obtain, after some manipulations taking also advantage of (3.8),

It is now more difficult to extract from (3.8)–(3.9) analytical information on the large time trends of the model and its equilibria due to the substantial lack of basic conservation properties. In section 4 we will explore numerically some representative case studies. Furthermore, we mention that in ^[24] a qualitative analysis is proposed concerning the aggregate time asymptotic trends of the quarantine model without network (i.e. formally $\left\lvert{ \mathcal{I}}\right\rvert = 1$) for both constant and non-constant label switch probabilities.

4.   Numerical tests

We present now some illustrative numerical tests, which exemplify the types of results and considerations that can be drawn from the models introduced in the previous sections. In each test, we compute the aggregate densities of individuals and the mean viral loads by solving both the stochastic particle model via suitable adaptations of classical Monte Carlo schemes for kinetic equations, cf. e.g., ^[25], and the macroscopic equations via standard numerical methods for ODEs such as the fourth-order Runge-Kutta method. More specifically, the particle models that we consider are (2.5), (2.22), (3.2), which are discrete in time and may therefore be straightforwardly translated into algorithmic schemes. For an explicit description of a Monte Carlo algorithm incorporating both particle interactions and label switching we refer to ^[24]. Compliance of the particle and macroscopic solutions validates the macroscopic models that we have obtained from the kinetic description of the particle models.

In all tests we consider the strongly connected graph illustrated in Figure 2, which models a network of e.g., three different cities. Hence $\mathcal{I} = \{1, \, 2, \, 3\}$. We choose the transition matrix as

and the parameters of the interactions as indicated in Table 1.

4.1.   Test 1: Basic interaction dynamics on a network

We begin from the basic model described in section 2.1, in which individuals move across the nodes of the network according to the transition probabilities encoded in the matrix (4.1) and exchange viral loads according to the interaction rules (2.3)–(2.4). The particle model is (2.5) while the corresponding macroscopic model is (2.9)–(2.12).

We choose the initial conditions as follows: in nodes $1$ and $2$ we set

while in node $3$ we let $f_{3, 0}(v)$ be an inverse-gamma distribution with shape parameter $3$, scale parameter $2$ and mass equal to $0.1$:

Therefore initially nodes $1$ and $2$ are disease-free, i.e. $m_{1, 0} = m_{2, 0} = 0$, while node $3$ features an infection onset with $m_{3, 0} = 1$. The initial densities are $\rho_{1, 0} = 0.3$, $\rho_{2, 0} = 0.6$, $\rho_{3, 0} = 0.1$.

Figure 3a, b show that the densities in the three nodes reach non-zero asymptotic values as predicted by Proposition 2.1 and that, since $\nu_1 < \nu_2$, the infection blows up in all nodes consistently with the analysis performed in section 2.2. These figures also show a perfect correspondence between the Monte Carlo solution of the particle model (2.5) and the macroscopic model (2.9)–(2.12).

Figure 3c, d show instead the effect of a cordon sanitaire around node $3$, which from time $t = 2$ onwards isolates that node from the others. We simulate this by modifying the transition matrix (4.1) for $t\geq 2$ as

In particular, $P_{33} = 1$ and $P_{i3} = P_{3j} = 0$ for $i, \, j = 1, \, 2$ because the cordon sanitaire prevents individuals from entering and exiting from node $3$. The densities reach new non-zero asymptotic values and the infection still blows up in all nodes because $\nu_1 < \nu_2$. Hence isolating the hotspot of the infection is per se useless as a confinement measure if it is not associated to other targeted interventions within the nodes. We will come back to this issue in section 4.3. For the moment, we observe however that the model predicts a different rate of blow-up of the infection in node $3$ compared to the case without cordon sanitaire. In particular, despite the fact that the infection originated from node $3$, the latter features the lowest blow-up rate, because it has been isolated from node $2$, which is the most populated and connected one. Again, we notice a nice matching between the solutions of the particle and the macroscopic models.

4.2.   Test 2: Commuters

We consider now the model of the commuters introduced in section 2.3, cf. also Figure 4. We choose the following commuting probabilities:

and we set $p_{ii} = 1$ for all $i = 1, \, 2, \, 3$. We recall that, as explained in Remark 2.2, the conceptual meaning of these probabilities is different from that of the transition probabilities (4.1). In particular, the latter are not used in this model.

At $t = 0$ we prescribe the following distributions:

meaning that initially all commuting routes are disease-free but the routes $3\to 2$ and $3\to 3$, where $m_{32, 0} = m_{33, 0} = 1$. Therefore node $3$ is again the infection hotspot. As a matter of fact, $3\to 3$ is actually not a real commuting route but identifies individuals who remain always in node $3$, i.e. do not commute.

The total masses of commuters initially present in the nodes are

while those of the commuters along the various routes are

As explained in section 2.3, these latter values are conserved in time. We also compute

noticing that $\bar{\rho}_{\{i, \, i\}} = 2\rho_{ii, 0}$ for all $i = 1, \, 2, \, 3$ for consistency with the definition of $\bar{\rho}_{\{i, \, j\}}$ given in section 2.3.

The particle model is now (2.22) and the macroscopic model is (2.24)–(2.26).

From Figure 5a we see that the total densities in the nodes quickly reach non-zero asymptotic values, which, plugging the numerical data above into (2.25), can be computed precisely as

We remark that, due to the repetitiveness of the commuting dynamics, it makes perfectly sense that the densities in the nodes are essentially constant in time apart from a short initial transient. We observe moreover that the gap $\left\lvert{\rho_i^\infty-\rho_{i, 0}}\right\rvert$ between the initial and equilibrium values is of order $10^{-2}$, hence reasonably small, for all $i = 1, \, 2, \, 3$.

From Figure 5b we observe once again a blow-up of the infection in all nodes like in Test 1. Nevertheless, unlike the case illustrated in Figure 3b, here the trends of the mean viral loads in the three nodes differ more consistently from each other. In particular, the blow-up occurs earlier in node $2$, i.e. the one more crossed by commuting routes, and later in node $1$, which is not only one of the two nodes less crossed by commuting routes but also an initially disease-free node. This shows that if one assumes more specific migration paths than simple random transitions from node to node then our models can provide more accurate and realistic predictions.

Finally, from Figure 5 we still observe a nice agreement between the microscopic particle dynamics and the aggregate macroscopic trends derived through the kinetic description.

4.3.   Test 3: Quarantine in the nodes

Finally we consider the model introduced in section 3 to explore the effect of quarantine as a measure to confine the infection node-wise. Remarkably, quarantine is one of the few confinement measures immediately implementable in case of new infectious diseases for which medical treatments are not yet available.

The network for this test is the one in Figure 3 with transition probabilities $P^1_{ij}$ of non-quarantined individuals coinciding with the corresponding $P_{ij}$'s in (4.1). The transition probabilities of quarantined individuals are instead the trivial ones $P^2_{ij} = \delta_{ij}$, consistently with the general setting presented in section 3. The other relevant parameters of the interactions are given in Table 1. In addition to them, we need to prescribe the label switch probabilities $T^{21}_i(v)$, $T^{12}_i(v)$, which describe the processes by which an individual with viral load $v$ in node $i$ is diagnosed as infected, hence put in isolation, and released from quarantine, hence readmitted in the society, respectively. We observe that $T^{21}_i(v)$, $T^{12}_i(v)$ may be correlated with the sensitivity of either testing technique used to detect the infection in an individual.

The initial conditions are the following:

They model a scenario in which at $t = 0$ no individual is quarantined in any node and node $3$ is the hotspot of the infection. In particular, we have

and

Furthermore, $m_{3, 0} = m^1_{3, 0} = 0.5$ while all other mean viral loads are initially zero.

We begin with the case of constant label switch probabilities: then the particle model (3.2) admits a macroscopic counterpart in closed form given by model (3.8)–(3.9). In particular, we let

Figure 6 shows the time trends of the total density and mean viral load (top row), as well as those of the density and mean viral load of non-quarantined and quarantined individuals (middle and bottom rows, respectively), in each node of the network. Apart from remarking again a perfect correspondence between the solutions to the particle and the macroscopic models, we observe from Figure 6b (and, with further specificity, from Figure 6d, f) that in the long run the infection is eradicated in every node because all mean viral loads tend to zero for $t\to +\infty$. It is interesting to compare this result with the one obtained in ^[24], where the same model with quarantine is analysed in the absence of a network. In practice, the model in ^[24] can be seen as the analogous of the present model in an isolated node, for instance one around which a cordon sanitaire has been established. Remarkably, in that case and with the very same values of the parameters used here an infection blow-up occurs (cf., [24,Figure 2]). Only if the probability for an individual to be diagnosed as infected and quarantined is sufficient higher than $0.2$ (the analysis in [24,sections 5.1,6] indicates that the minimal threshold is $0.28$) the infection may be successfully kept under control and eradicated in the long run. The comparison between these two results puts in evidence a striking effect of the network: allowing the individuals to migrate and mix in different locations, with all other features unchanged including the effectiveness of the testing techniques, may help dominate the infection more effectively than the confinement produced by a cordon sanitaire. This conclusion, here emerging naturally as a consequence of our quite general model, is in agreement with the observations made in ^[27] using more ad-hoc compartmental models specifically conceived to address this phenomenon.

To conclude, we further consider the case of non-constant label switch probabilities. Specifically, we set

Hence the probability for an individual to be diagnosed as infected and quarantined increases monotonically when the viral load increases, while the probability to be recognised as recovered and readmitted in the society decreases monotonically. This is more consistent than constant label switch probabilities with the way in which real screening tests for infection detection work, see e.g., ^[19]. Notice that these functions are bounded from above by the values formerly used as constant label switch probabilities:

For this case we have not derived a macroscopic model in closed form (though some hints in this direction may be found in ^[24] in the networkless setting), therefore in Figure 7 we only show the results of the particle model (3.2). We notice in particular that the quarantine is still successful in keeping the infection globally under control, indeed no blow-up of the mean viral load occurs in any node. Nevertheless, since the label switch probabilities are not the same for all individuals, unlike the previous case the infection cannot be eradicated by relying only on quarantine. The infection becomes instead endemic as shown by the fact that the mean viral loads reach asymptotically comparable non-zero values in all nodes, cf., Figure 7b. The details of Figure 7d, f suggest furthermore that the mean viral loads of the non-quarantined individuals set reasonably on lower endemic values than those of the quarantined individuals. On the other hand, Figure 7c, e show that in the endemic regime the number of non-quarantined individuals is systematically higher than that of quarantined ones in all nodes, a fact that meets again the intuitive expectation. The lower number $\rho^2_i$ of quarantined individuals in all nodes also explains why the trend of the global mean viral loads $m_i$ (Figure 7b) is quite close to that of the mean viral loads of non-quarantined individuals $m^1_i$ (Figure 7d) for all $i = 1, \, 2, \, 3$, indeed we recall that $m_i = \frac{\rho^1_i}{\rho_i}m^1_i+\frac{\rho^2_i}{\rho_i}m^2_i$.

By suitably increasing (4.2), i.e. the sensitivity of the testing techniques, one may expect that the model predicts the eradication of the infection in the long run also in this case.

5.   Conclusions

In this paper, we have proposed a formal derivation of population models which describe the spreading of an infectious disease on a spatial network by taking into account the role of the viral load of the individuals. In particular, with the introduction of the viral load as a descriptive variable of the system we have modelled the transmission of the disease in a more specific way than by simply estimating the number of infectious contacts out of the gross number of susceptible and infectious individuals. Furthermore, we have avoided the necessity to partition the population of each node in infection-dependent compartments, because we have described both the individual and the aggregate epidemiological trends by means of the individual and mean viral loads, respectively.

Our derivation has taken advantage of the conceptual paradigms of statistical mechanics and kinetic theory. Starting from a particle model, in which individuals are characterised by a microscopic state comprising their viral load and their current node in the network, we have provided a mesoscopic description in terms of Boltzmann-type equations for the distributions of the viral load in the various nodes. Subsequently, thanks to the linear structure of the interaction rules modelling the pairwise transmission of viral load and the jumps across the nodes of the network, we have been able to obtain closed systems of non-linear macroscopic equations for the time evolution of the density of individuals and their mean viral load in each node. We have then characterised analytically the large-time aggregate trends, such as the existence of equilibrium density distributions and the blow-up or eradication of the infection in the nodes, in some representative cases in terms of relevant microscopic parameters. Finally, we have also extended the basic model of disease transmission on networks to even more realistic scenarios, such as the one of commuters, who move across the nodes according to more specific criteria than simply random jumps, and the one in which quarantine is applied in the nodes as a confinement measure of the infection. In all these cases we have obtained the macroscopic counterpart of the particle model via a statistical mechanics and kinetic theory approach. Our numerical tests have confirmed the matching between the particle and the macroscopic models, thereby validating the latter as reliable approximations of the former more amenable to analytical investigations and quick and accurate numerical solutions. As a by-product, our approach has provided a contribution to the broader topic of kinetic models of network-structured interactions by addressing Boltzmann-type kinetic equations on graphs.

Deliberately, we have not tried to match real scenarios observed during a pandemic by calibrating or comparing the results of our models with data. The present work is indeed a methodological one, specifically focused on a rigorous formal derivation of new population models, which can describe aspects normally neglected in standard epidemiological models. Accordingly, all parameters and variables of our models have to be understood as dimensionless. Indeed, in our analyses the main point is not their numerical values but rather the relationships among them and their orders of magnitude, which may induce different regimes and different asymptotic scenarios. In this respect, our tests have essentially explored plausible prototypical situations, while the systematic application of our models to real case studies will be the object of future investigations. We remark however that to adapt our models to the study of real-life cases it is in principle sufficient to identify dimensional reference values for each parameter and variable and scale the non-dimensional quantities accordingly. We also point out that, although in a qualitative dimensionless form, our results already indicate the ability of our models to address interesting issues typically out of the scope of standard compartmental models. On one hand, the presence of the spatial network has reproduced the observed possible inefficiency of mobility restrictions to control the growth of epidemics, cf., ^[27]. On the other hand, the introduction of the viral load in the microscopic state of the individuals makes it possible to investigate quantitatively the interplay between sensitivity and frequency of the viral tests for an optimal screening of the population. As witnessed by the contemporary literature, see e.g., ^[19], this is a particularly relevant issue in the fight against newly discovered infectious disease, however still addressed in a largely qualitative and empirical way.

Acknowledgements

This research was partially supported by the Italian Ministry for Education, University and Research (MIUR) through the "Dipartimenti di Eccellenza" Programme (2018-2022), Department of Mathematical Sciences "G. L. Lagrange", Politecnico di Torino (CUP: E11G18000350001) and through the PRIN 2017 project (No. 2017KKJP4X) "Innovative numerical methods for evolutionary partial differential equations and applications".

NL is a postdoctoral research fellow ("titolare di Assegno di Ricerca") of Istituto Nazionale di Alta Matematica (INdAM, Italy).

Both authors are members of GNFM (Gruppo Nazionale per la Fisica Matematica) of INdAM, Italy.