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

Nadia Loy<EMAIL_ADDRESS>Andrea Tosin<EMAIL_ADDRESS>
(Department of Mathematical Sciences “G. L. Lagrange”
Politecnico di Torino, Italy)

###### Abstract

In this paper, we propose a Boltzmann-type kinetic model of the spreading of
an infectious disease on a network. The latter describes the connections among
countries, cities or districts depending on the spatial scale of interest. The
disease transmission is represented in terms of the viral load of the
individuals and is mediated by social contacts among them, taking into account
their displacements across the nodes of the network. We formally derive the
hydrodynamic equations for the density and the mean viral load of the
individuals on the network and we analyse the large-time trends of these
quantities with special emphasis on the cases of blow-up or eradication of the
infection. By means of numerical tests, we also investigate the impact of
confinement measures, such as quarantine or localised lockdown, on the
diffusion of the disease on the network.

Keywords: Boltzmann-type equations, Markov-type jump processes, label
switching, commuters, quarantine

Mathematics Subject Classification: 35Q20, 35Q70, 35Q84

##  1 Introduction

The spreading of infectious diseases is intimately correlated with patterns of
human movement, as clearly demonstrated by empirical observations [22]. Yet,
many compartmental epidemiological models, inspired by the celebrated Kermack-
McKendrick theory [15], 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., [8, 14]. 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., [1, 6] and [18, 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 may move from one node to another the description is on the whole
spatially inhomogeneous and appears to be particularly suited to model main
and preferential paths. 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. [2, 3, 4, 8]. 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 recent survey [23] and the
specific application [21] 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 [12, 13].

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., [9,
10], 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
[16]. 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 models of social interactions, see e.g., [7] for a very
recent proposal. 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
[17] 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

Figure 1: A strongly connected graph in which individuals interact within the
vertices and may jump from vertex to vertex.

###  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

$P_{ij}:=\operatorname{Prob}(j\to i)\in[0,\,1],\qquad i,\,j\in\mathcal{I},$
(1)

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

$\sum_{i\in\mathcal{I}}P_{ij}=1,\qquad\forall\,j\in\mathcal{I},$ (2)

hence the square matrix
$\mathbf{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 $\mathbf{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 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

$v^{\prime}=pv+qv_{\ast},\qquad v_{\ast}^{\prime}=pv_{\ast}+qv,$ (3)

where $v,\,v_{\ast}$ represent the pre-interaction viral loads of two
interacting individuals and $v^{\prime},\,v_{\ast}^{\prime}$ their post-
interaction viral loads. Moreover, $p,\,q\in\mathbb{R}_{+}$ are either
deterministic or stochastic interaction coefficients. The binary rules (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 (3), whereas
individuals in different vertices do not interact. The reason is that
contagion dynamics require a certain proximity to be effective. Notice that
(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

$p=1-\nu_{1}+\eta,\qquad q=\nu_{2},$ (4)

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^{\prime},\,v_{\ast}^{\prime}\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:

$\displaystyle\begin{aligned} X_{t+\Delta{t}}&=(1-\Theta)X_{t}+\Theta
J_{t},\\\
V_{t+\Delta{t}}&=\left(1-\Xi\delta_{X_{t},X^{\ast}_{t}}\right)V_{t}+\Xi\delta_{X_{t},X^{\ast}_{t}}V^{\prime}_{t},\end{aligned}$
(5)

where:

  1. 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

$\operatorname{Prob}(\Theta=1)=\chi\Delta{t},\qquad\operatorname{Prob}(\Xi=1)=\mu\Delta{t},$
(6)

$\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;

  2. ii)

$J_{t}\in\mathcal{I}$ is a random variable returning the new vertex after a
vertex jump, with

$\operatorname{Prob}(J_{t}=i|X_{t}=j)=P_{ij},\qquad i,\,j\in\mathcal{I};$

  3. iii)

$V^{\prime}_{t}\in\mathbb{R}_{+}$ is the new viral load after a social
contact. In view of (3) we have $V^{\prime}_{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;

  4. iv)

$\delta_{X_{t},X^{\ast}_{t}}$ is the Kronecker delta:

$\delta_{X_{t},X^{\ast}_{t}}:=\begin{cases}1&\text{if }X_{t}=X^{\ast}_{t}\\\
0&\text{otherwise},\end{cases}$

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 (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 [17], see also [20] 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:

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f_{i}(v,t)\,dv&=\chi\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{j\in\mathcal{I}}P_{ij}f_{j}(v,t)-f_{i}(v,t)\right)dv\\\
&\phantom{=}+\mu\int_{\mathbb{R}_{+}}\varphi(v)Q(f_{i})(v,t)\,dv,&i\in\mathcal{I},\end{aligned}$
(7)

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:

$\int_{\mathbb{R}_{+}}\varphi(v)Q(f_{i})(v,t)\,dv=\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}}\langle\varphi(v^{\prime})-\varphi(v)\rangle
f_{i}(v,t)f_{i}(v_{\ast},t)\,dv\,dv_{\ast},$ (8)

where $v^{\prime}$ is given by (3)-(4) as a function of the integration
variables $v$, $v_{\ast}$.

###  2.2 Aggregate description of mobility and contagion

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

$\rho_{i}(t):=\int_{\mathbb{R}_{+}}f_{i}(v,t)\,dv\geq 0,$

is in general not conserved in time because of the jumps across the vertices.
Setting $\varphi(v)=1$ in (7) we obtain that $\rho_{i}$ satisfies the equation

$\frac{d\rho_{i}}{dt}=\chi\left(\sum_{j\in\mathcal{I}}P_{ij}\rho_{j}-\rho_{i}\right),\qquad
i\in\mathcal{I},$ (9)

which in vector notation reads

$\frac{d\boldsymbol{\rho}}{dt}=\chi(\mathbf{P}-\mathbf{I})\boldsymbol{\rho}$
(10)

with $\boldsymbol{\rho}:=(\rho_{i})_{i\in\mathcal{I}}$. From (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

$(\mathbf{P}-\mathbf{I})\boldsymbol{\rho}^{\infty}=\mathbf{0}.$ (11)

In practice, $\boldsymbol{\rho}^{\infty}$ should be an eigenvector of
$\mathbf{P}$ corresponding to the eigenvalue $1$. Notice that $\mathbf{P}$
admits indeed this eigenvalue, which coincides also with its spectral radius,
because it is a stochastic matrix. Moreover, since $\mathbf{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., [19] 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 (9) and (2) we deduce

$\frac{d}{dt}\|\boldsymbol{\rho}(t)\|_{1}:=\frac{d}{dt}\sum_{i\in\mathcal{I}}\rho_{i}(t)=0.$

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 (11).

Still from Perron-Frobenius theory we have that $1$ is also the spectral limit
of $\mathbf{P}$, namely the maximum real part of its eigenvalues. Therefore
the real parts of the eigenvalues of the matrix $\mathbf{P}-\mathbf{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 (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 (11), which is a stable
and attractive asymptotic density distribution for (10).

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

$m_{i}(t):=\frac{1}{\rho_{i}(t)}\int_{\mathbb{R}_{+}}vf_{i}(v,t)\,dv,$

which, recalling also (3)-(4), turns out to satisfy the equation

$\frac{d}{dt}(\rho_{i}m_{i})=\chi\left(\sum_{j\in\mathcal{I}}P_{ij}\rho_{j}m_{j}-\rho_{i}m_{i}\right)+\mu(\nu_{2}-\nu_{1})\rho_{i}^{2}m_{i}.$
(12)

From here, invoking (2) we discover

$\frac{d}{dt}\sum_{i\in\mathcal{I}}\rho_{i}m_{i}=\mu(\nu_{2}-\nu_{1})\sum_{i\in\mathcal{I}}\rho_{i}^{2}m_{i},$

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
$\mathbf{m}^{\infty}=\mathbf{0}$, which represents the eradication of the
infection in all nodes of the network. Linearising the previous equation
around the equilibrium configuration
$(\boldsymbol{\rho},\,\mathbf{m})=(\boldsymbol{\rho}^{\infty},\,\mathbf{0})$
we obtain

$\frac{d}{dt}\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i}=\mu(\nu_{2}-\nu_{1})\sum_{i\in\mathcal{I}}(\rho_{i}^{\infty})^{2}m_{i}.$
(13)

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

$\sum_{i\in\mathcal{I}}(\rho_{i}^{\infty})^{2}m_{i}\geq
c\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i}$

to conclude that:

  1. i)

if $\nu_{1}>\nu_{2}$ then

$\frac{d}{dt}\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i}\leq
c\mu(\nu_{2}-\nu_{1})\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i},$

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
$\mathbf{m}^{\infty}=\mathbf{0}$ is locally asymptotically stable, hence the
infection may be globally eradicated in the long run;

  2. ii)

if $\nu_{1}<\nu_{2}$ then

$\frac{d}{dt}\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i}\geq
c\mu(\nu_{2}-\nu_{1})\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i},$

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
$\mathbf{m}^{\infty}=\mathbf{0}$ is unstable: the infection is not under
control and tends to grow in time;

  3. iii)

if $\nu_{1}=\nu_{2}$ then $\sum_{i\in\mathcal{I}}\rho_{i}m_{i}$ is constant in
time. This indicates that the configuration $\mathbf{m}^{\infty}=\mathbf{0}$
is stable but not attractive: the infection is under control but cannot be
eradicated and becomes endemic.

####  2.2.1 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 (13) modifies as

$\frac{d}{dt}\sum_{i\in\mathcal{I}}\rho_{i}^{\infty}m_{i}=\mu\sum_{i\in\mathcal{I}}(\nu_{2,i}-\nu_{1,i})(\rho_{i}^{\infty})^{2}m_{i}.$

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 (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

$\chi=1,\qquad\mu=\frac{1}{\epsilon},$

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 (7) is possible in:

  1. i)

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

$\frac{d}{d\tau}\int_{\mathbb{R}_{+}}\varphi(v)f_{i}(v,\tau)\,dv=\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}}\langle\varphi(v^{\prime})-\varphi(v)\rangle
f_{i}(v,\tau)f_{i}(v_{\ast},\tau)\,dv\,dv_{\ast},$ (14)

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;

  2. ii)

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

$\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)M_{i}(v,t)\,dv=\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{j\in\mathcal{I}}P_{ij}M_{j}(v,t)-M_{i}(v,t)\right)dv,$
(15)

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

Taking $\varphi(v)=1$ in (14) we discover that in all vertices the density
$\rho_{i}$ is constant on the time scale $\tau$. Taking instead $\varphi(v)=v$
we get

$\frac{dm_{i}}{d\tau}=(\nu_{2,i}-\nu_{1,i})\rho_{i}m_{i}.$ (16)

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 (16) 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
(15) and choosing $\varphi(v)=1$ yields the same time evolution of the
densities $\rho_{i}$, $i\in\mathcal{I}$, as (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:

$\frac{d}{dt}(\rho_{i^{\ast}}m_{i^{\ast}})=(P_{i^{\ast}i^{\ast}}-1)\rho_{i^{\ast}}m_{i^{\ast}}.$

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 (15) with
$\varphi(v)=v$ we deduce

$\frac{d}{dt}(\rho_{i}m_{i})=\sum_{j\in\mathcal{I}\setminus\\{i^{\ast}\\}}P_{ij}\rho_{j}m_{j}-\rho_{i}m_{i},\qquad
i\in\mathcal{I}\setminus\\{i^{\ast}\\},$

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

$\displaystyle\begin{aligned}
\frac{d}{dt}\sum_{i\in\mathcal{I}\setminus\\{i^{\ast}\\}}\rho_{i}m_{i}&=\sum_{j\in\mathcal{I}\setminus\\{i^{\ast}\\}}\left(\sum_{i\in\mathcal{I}\setminus\\{i^{\ast}\\}}P_{ij}\right)\rho_{j}m_{j}-\sum_{i\in\mathcal{I}\setminus\\{i^{\ast}\\}}\rho_{i}m_{i}\\\
&=-\sum_{j\in\mathcal{I}\setminus\\{i^{\ast}\\}}P_{i^{\ast}j}\rho_{j}m_{j},\end{aligned}$
(17)

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:

$\mathcal{I}^{\ast}:=\\{j\in\mathcal{I}\setminus\\{i^{\ast}\\}\,:\,P_{i^{\ast}j}>0\\}.$

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 (17) we deduce

$\frac{d}{dt}\sum_{i\in\mathcal{I}\setminus\\{i^{\ast}\\}}\rho_{i}m_{i}\leq-a\sum_{j\in\mathcal{I}^{\ast}}\rho_{j}m_{j},$

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

$\sum_{i\in\mathcal{I}^{\ast}}\rho_{i}(t)m_{i}(t)\leq\sum_{i\in\mathcal{I}\setminus\\{i^{\ast}\\}}\rho_{i,0}m_{i,0}-a\int_{0}^{t}\sum_{j\in\mathcal{I}^{\ast}}\rho_{j}(s)m_{j}(s)\,ds.$

By Grönwall’s inequality we deduce then

$\sum_{i\in\mathcal{I}^{\ast}}\rho_{i}(t)m_{i}(t)\leq
e^{-at}\sum_{i\in\mathcal{I}\setminus\\{i^{\ast}\\}}\rho_{i,0}m_{i,0},$

which says that $\sum_{i\in\mathcal{I}^{\ast}}\rho_{i}m_{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 (3) with coefficients given by (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 (16) 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

$\displaystyle\begin{aligned}
(X_{o,t+\Delta{t}},\,X_{d,t+\Delta{t}})&=(1-\Theta)(X_{o,t},\,X_{d,t})+\Theta(X^{\prime}_{o,t},\,X^{\prime}_{d,t})\\\
V_{t+\Delta{t}}&=\left(1-\Xi\delta_{X_{o,t},X^{\ast}_{o,t}}\right)V_{t}+\Xi\delta_{X_{o,t},X^{\ast}_{o,t}}V^{\prime}_{t},\end{aligned}$
(18)

where $\Theta$, $\Xi$ are like in (6), $V^{\prime}_{t}=pV_{t}+qV^{\ast}_{t}$
from (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^{\prime}_{o,t},\,X^{\prime}_{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^{\prime}_{o,t}=X_{d,t}$ and
$X^{\prime}_{d,t}=X_{o,t}$ if the journey takes place or
$X^{\prime}_{o,t}=X_{o,t}$ and $X^{\prime}_{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^{\prime}_{o,t},\,X^{\prime}_{d,t})$ is

$\displaystyle\operatorname{Prob}((X^{\prime}_{o,t},\,X^{\prime}_{d,t})=(X_{d,t},\,X_{o,t}))$
$\displaystyle=p_{X_{d,t},X_{o,t}},$
$\displaystyle\operatorname{Prob}((X^{\prime}_{o,t},\,X^{\prime}_{d,t})=(X_{o,t},\,X_{d,t}))$
$\displaystyle=1-p_{X_{d,t},X_{o,t}}.$

We are now in a position to deduce from (18) 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 [17] with
minimal adaptations, in the continuous-time limit $\Delta{t}\to 0^{+}$ we
formally obtain from (18) that $f_{ji}$ satisfies the following Boltzmann-type
kinetic equation in weak form:

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f_{ji}(v,t)\,dv&=\chi\int_{\mathbb{R}_{+}}\varphi(v)\left(p_{ij}f_{ij}(v,t)-p_{ji}f_{ji}(v,t)\right)dv\\\
&\phantom{=}+\mu\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}}\langle\varphi(v^{\prime})-\varphi(v)\rangle
f_{ji}(v,t)f_{i}(v_{\ast},t)\,dv\,dv_{\ast},&i,\,j\in\mathcal{I},\end{aligned}$
(19)

where $v^{\prime}$ is given by (3)-(4) and

$f_{i}(v,t):=\sum_{j\in\mathcal{I}}f_{ji}(v,t)$

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 (19) features two important differences with respect to (7), that we
now elucidate.

  1. i)

The commuting term

$\int_{\mathbb{R}_{+}}\varphi(v)\left(p_{ij}f_{ij}(v,t)-p_{ji}f_{ji}(v,t)\right)dv$

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}$.

  2. ii)

The contagion term

$\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}}\langle\varphi(v^{\prime})-\varphi(v)\rangle
f_{ji}(v,t)f_{i}(v_{\ast},t)\,dv\,dv_{\ast}$

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. (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

$\rho_{ji}(t):=\int_{\mathbb{R}_{+}}f_{ji}(v,t)\,dv$

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

$\frac{d\rho_{ji}}{dt}=\chi\left(p_{ij}\rho_{ij}-p_{ji}\rho_{ji}\right),\qquad
i,\,j\in\mathcal{I}.$ (20)

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 by111We 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.
$\bar{\rho}_{\\{i,j\\}}:=\rho_{ji}+\rho_{ij}$ the total mass of
$ij$-commuters. Then $\rho_{ij}=\bar{\rho}_{\\{i,j\\}}-\rho_{ji}$, whence

$\frac{d\rho_{ji}}{dt}=\chi\left(p_{ij}\bar{\rho}_{\\{i,j\\}}-(p_{ji}+p_{ij})\rho_{ji}\right),$

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

$\rho^{\infty}_{ji}:=\frac{p_{ij}}{p_{ji}+p_{ij}}\bar{\rho}_{\\{i,j\\}},\qquad\rho^{\infty}_{ij}:=\frac{p_{ji}}{p_{ji}+p_{ij}}\bar{\rho}_{\\{i,j\\}},\qquad
i,\,j\in\mathcal{I}.$

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

$\rho_{i}(t):=\sum_{j\in\mathcal{I}}\rho_{ji}(t)=\int_{\mathbb{R}_{+}}f_{i}(v,t)\,dv,$

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

$\rho^{\infty}_{i}=\sum_{j\in\mathcal{I}}\frac{p_{ij}}{p_{ji}+p_{ij}}\bar{\rho}_{\\{i,j\\}},\qquad
i\in\mathcal{I}.$ (21)

######  Remark 2.3.

From (20), the evolution equation for the total density $\rho_{i}$ turns out
to be

$\frac{d\rho_{i}}{dt}=\chi\sum_{j\in\mathcal{I}}\left(p_{ij}\rho_{ij}-p_{ji}\rho_{ji}\right),$

which is here the counterpart of (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

$\frac{d\rho_{i}}{dt}=\chi\left(\sum_{j\in\mathcal{I}}\rho_{ij}-\rho_{i}\right),\qquad
i\in\mathcal{I},$

which reminds more closely of (9) but with $P_{ij}\rho_{j}$ replaced by
$\rho_{ij}$. Notice that, unlike (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

$m_{ji}(t):=\frac{1}{\rho_{ji}(t)}\int_{\mathbb{R}_{+}}vf_{ji}(v,t)\,dv$

the mean viral load of $ij$-commuters in vertex $i$ at time $t$ then from (19)
with $\varphi(v)=v$ and taking also (3)-(4) and (20) into account we obtain

$\frac{dm_{ji}}{dt}=\chi
p_{ij}\frac{\rho_{ij}}{\rho_{ji}}(m_{ij}-m_{ji})+\mu\rho_{i}\left(\nu_{2}m_{i}-\nu_{1}m_{ji}\right),\qquad
i,\,j\in\mathcal{I},$ (22)

where

$m_{i}(t):=\frac{1}{\rho_{i}(t)}\int_{\mathbb{R}_{+}}vf_{i}(v,t)\,dv=\frac{1}{\rho_{i}(t)}\sum_{j\in\mathcal{I}}\rho_{ji}(t)m_{ji}(t)$

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 (22)).
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
[17], 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

$T^{hk}_{i}(v):=\operatorname{Prob}(k\to h|i,\,v)\in[0,\,1],\qquad
h,\,k\in\\{1,\,2\\},\ i\in\mathcal{I},\ v\in\mathbb{R}_{+},$

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

$T^{1k}_{i}(v)+T^{2k}_{i}(v)=1\qquad\forall\,k\in\\{1,\,2\\},\
i\in\mathcal{I},\ v\in\mathbb{R}_{+}.$ (23)

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

$\displaystyle\begin{aligned} X_{t+\Delta{t}}&=(1-\Theta)X_{t}+\Theta J_{t}\\\
H_{t+\Delta{t}}&=(1-\Lambda)H_{t}+\Lambda I_{t}\\\
V_{t+\Delta{t}}&=\left(1-\Xi\right)V_{t}+\Xi V^{\prime}_{t},\end{aligned}$
(24)

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

$\operatorname{Prob}(\Lambda=1)=\lambda\Delta{t}$

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

$\operatorname{Prob}(I_{t}=h|H_{t}=k,\,X_{t}=i,\,V_{t}=v)=T^{hk}_{i}(v).$

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^{\prime}_{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 (3). To take these facts into account we modify the law of $J_{t}$ as
follows:

$\operatorname{Prob}(J_{t}=i|X_{t}=j,\,H_{t}=h)=P^{h}_{ij}\in[0,\,1],$

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

$P^{2}_{ij}=\delta_{ij},\qquad\forall\ i,\,j\in\mathcal{I}$

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

$v^{\prime}=(1-\gamma+\xi)v,$ (25)

which is of the form (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^{\prime}\geq 0$
for all $v\geq 0$. On the whole, we set

$V^{\prime}_{t}=\begin{cases}(1-\nu_{1}+\eta)V_{t}+\nu_{2}V^{\ast}_{t}&\text{if
}H_{t}=H^{\ast}_{t}=1\text{ and }X_{t}=X^{\ast}_{t}\\\ V_{t}&\text{if
}(H_{t}=H^{\ast}_{t}=1\text{ and }X_{t}\neq X^{\ast}_{t})\text{ or
}(H_{t}=1\text{ and }H^{\ast}_{t}=2)\\\ (1-\gamma+\xi)V_{t}&\text{if
}H_{t}=2,\end{cases}$

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 (7) and the results in [17] we deduce from (24)
that in the continuous-time limit $\Delta{t}\to 0^{+}$ the $f^{h}_{i}$’s
satisfy formally the following kinetic equation in weak form:

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f^{h}_{i}(v,t)\,dv&=\chi\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{j\in\mathcal{I}}P^{h}_{ij}f^{h}_{j}(v,t)-f^{h}_{i}(v,t)\right)dv\\\
&\phantom{=}+\lambda\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{k=1}^{2}T^{hk}_{i}(v)f^{k}_{i}(v,t)-f^{h}_{i}(v,t)\right)dv\\\
&\phantom{=}+\mu\int_{\mathbb{R}_{+}}\varphi(v)Q^{h}(f^{h}_{i})(v,t)\,dv,&i\in\mathcal{I},\,h\in\\{1,\,2\\},\end{aligned}$
(26)

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 (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$):

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f^{1}_{i}(v,t)\,dv&=\chi\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{j\in\mathcal{I}}P^{1}_{ij}f^{1}_{j}(v,t)-f^{1}_{i}(v,t)\right)dv\\\
&\phantom{=}+\lambda\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{k=1}^{2}T^{1k}_{i}(v)f^{k}_{i}(v,t)-f^{1}_{i}(v,t)\right)dv\\\
&\phantom{=}+\mu\int_{\mathbb{R}_{+}}\int_{\mathbb{R}_{+}}\langle\varphi(v^{\prime})-\varphi(v)\rangle
f^{1}_{i}(v,t)f^{1}_{i}(v_{\ast},t)\,dv\,dv_{\ast},&i\in\mathcal{I}\end{aligned}$
(27)

with $v^{\prime}$ given by (3)-(4);

  * •

quarantined individuals ($h=2$):

$\displaystyle\begin{aligned}
\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(v)f^{2}_{i}(v,t)\,dv&=\lambda\int_{\mathbb{R}_{+}}\varphi(v)\left(\sum_{k=1}^{2}T^{2k}_{i}(v)f^{k}_{i}(v,t)-f^{2}_{i}(v,t)\right)dv\\\
&\phantom{=}+\mu\int_{\mathbb{R}_{+}}\langle\varphi(v^{\prime})-\varphi(v)\rangle
f^{2}_{i}(v,t)\,dv,&i\in\mathcal{I}\end{aligned}$ (28)

with $v^{\prime}$ given by (25).

###  3.1 Hydrodynamic model with quarantine on the network

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

$\rho^{h}_{i}(t):=\int_{\mathbb{R}_{+}}f^{h}_{i}(v,t)\,dv$

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

$\frac{d\rho^{h}_{i}}{dt}=\chi\left(\sum_{j\in\mathcal{I}}P^{h}_{ij}\rho^{h}_{j}-\rho^{h}_{i}\right)+\lambda\left(\sum_{k=1}^{2}\int_{\mathbb{R}_{+}}T^{hk}_{i}(v)f^{k}_{i}(v,t)\,dv-\rho^{h}_{i}\right),$
(29)

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

######  Remark 3.1.

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

$\rho_{i}(t):=\sum_{h=1}^{2}\rho^{h}_{i}(t)$

we find

$\frac{d\rho_{i}}{dt}=\chi\left(\sum_{j\in\mathcal{I}}\sum_{h=1}^{2}P^{h}_{ij}\rho^{h}_{j}-\rho_{i}\right),$

where we have recalled (23). If the transition probabilities from vertex to
vertex do not depend on $h$ then this equation coincides with (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 (29) 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, (29) becomes instead
a self-consistent equation for the densities $\rho^{h}_{i}$:

$\frac{d\rho^{h}_{i}}{dt}=\chi\left(\sum_{j\in\mathcal{I}}P^{h}_{ij}\rho^{h}_{j}-\rho_{i}\right)+\lambda\left(\sum_{k=1}^{2}T^{hk}_{i}\rho^{k}_{i}-\rho^{h}_{i}\right),$

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

$\begin{cases}\dfrac{d\rho^{1}_{i}}{dt}=\chi\left(\displaystyle{\sum_{j\in\mathcal{I}}}P^{1}_{ij}\rho^{1}_{j}-\rho^{1}_{i}\right)+\lambda\left(T^{12}_{i}\rho^{2}_{i}-T^{21}_{i}\rho^{1}_{i}\right)\\\\[14.22636pt]
\dfrac{d\rho^{2}_{i}}{dt}=\lambda\left(T^{21}_{i}\rho^{1}_{i}-T^{12}_{i}\rho^{2}_{i}\right),\end{cases}\quad
i\in\mathcal{I}.$ (30)

System (30) 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 (27)-(28) to deduce a coupled system of
hydrodynamic equations for the mean viral loads

$m^{h}_{i}(t):=\frac{1}{\rho^{h}_{i}(t)}\int_{\mathbb{R}_{+}}vf^{h}_{i}(v,t)\,dv,\qquad
i\in\mathcal{I},\,h\in\\{1,\,2\\}.$

Letting $\varphi(v)=v$ in (27)-(28) and assuming again that the $T^{hk}_{i}$’s
are constant with respect to $v$ we obtain, after some manipulations taking
also advantage of (30),

$\begin{cases}\dfrac{dm^{1}_{i}}{dt}=\chi\displaystyle{\sum_{j\in\mathcal{I}}}P^{1}_{ij}\dfrac{\rho^{1}_{j}}{\rho^{1}_{i}}\left(m^{1}_{j}-m^{1}_{i}\right)+\lambda
T^{12}_{i}\dfrac{\rho^{2}_{i}}{\rho^{1}_{i}}\left(m^{2}_{i}-m^{1}_{i}\right)+\mu(\nu_{2}-\nu_{1})\rho^{1}_{i}m^{1}_{i}\\\\[14.22636pt]
\dfrac{dm^{2}_{i}}{dt}=\lambda
T^{21}_{i}\dfrac{\rho^{1}_{i}}{\rho^{2}_{i}}\left(m^{1}_{i}-m^{2}_{i}\right)-\mu\gamma
m^{2}_{i},\end{cases}\quad i\in\mathcal{I}.$ (31)

It is now more difficult to extract from (30)-(31) 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
[17] 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

Table 1: Parameters of the interactions used in the numerical tests of Section 4. The symbol $\mathcal{U}([-a,\,a])$ indicates a random variable uniformly distributed in the interval $[-a,\,a]\subset\mathbb{R}$, $a>0$ Parameter | Test 1 | Test 2 | Test 3
---|---|---|---
(Figure 3) | (Figure 5) | (Figures 6, 7)
$\chi$ | $1$ | $1$ | $1$
$\mu$ | $1$ | $1$ | $1$
$\lambda$ | – | – | $1$
$\nu_{1}$ | $0$ | $0$ | $0$
$\nu_{2}$ | $0.2$ | $0.2$ | $0.2$
$\eta$ | $\mathcal{U}([-1,\,1])$ | $\mathcal{U}([-1,\,1])$ | $\mathcal{U}([-1,\,1])$
$\gamma$ | – | – | $0.3$
$\xi$ | – | – | $\mathcal{U}([-0.7,\,0.7])$
Figure 2: The three-node network considered in the numerical tests of Section
4

We present now some illustrative numerical tests, which exemplify the type 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 particle model via
suitable adaptations of classical Monte Carlo schemes for kinetic equations,
cf. e.g., [20], and the macroscopic equations via standard numerical methods
for ODEs such as the fourth-order Runge-Kutta method. Compliance of the two
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

$\mathbf{P}=\begin{pmatrix}P_{11}&P_{12}&P_{13}\\\ P_{21}&P_{22}&P_{23}\\\
P_{31}&P_{32}&P_{33}\end{pmatrix}=\begin{pmatrix}0.3&0.1&0\\\ 0.7&0.7&0.5\\\
0&0.2&0.5\end{pmatrix}$ (32)

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

###  4.1 Test 1: Basic interaction dynamics on a network

Figure 3: Test 1 (basic interaction dynamics on a network, Section 4.1). Top
row: standard evolution. Bottom row: cordon sanitaire around node $3$ for
$t\geq 2$

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 (32) and exchange viral loads according to the
interaction rules (3)-(4). The particle model is (5) while the corresponding
macroscopic model is (9)-(12).

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

$f_{1,0}(v)=0.3\delta(v),\qquad f_{2,0}(v)=0.6\delta(v)$

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$:

$f_{3,0}(v)=\frac{2e^{-\frac{2}{v}}}{5v^{4}}$

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$.

Figures 3a, 3b 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 (5) and
the macroscopic model (9)-(12).

Figures 3c, 3d 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 (32) for $t\geq 2$ as

$\mathbf{P}=\begin{pmatrix}0.3&0.3&0\\\ 0.7&0.7&0\\\ 0&0&1\end{pmatrix}.$

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 lower rate of blow-up of the infection in all nodes 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

Figure 4: The three-node network for the problem of the commuters discussed in
Section 4.2. Blue particles travel only between nodes $1$ and $2$ while red
particles travel only between nodes $2$ and $3$. Moreover, in every node there
is a portion of particles that do not commute

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

$\begin{array}[]{ll}p_{21}=0.7,&p_{31}=0,\\\ p_{12}=0.3,&p_{32}=0.8,\\\
p_{13}=0,&p_{23}=0.2,\end{array}$

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 (32). In particular, the latter are not
used in this model.

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

$\begin{array}[]{lll}f_{11,0}(v)=0.21\delta(v),&f_{12,0}(v)=0.06\delta(v),&f_{13,0}(v)=0,\\\
f_{21,0}(v)=0.09\delta(v),&f_{22,0}(v)=0.3\delta(v),&f_{23,0}(v)=0.05e^{-v},\\\
f_{31,0}(v)=0,&f_{32,0}(v)=0.24\delta(v),&f_{33,0}(v)=0.05e^{-v},\end{array}$

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

$\rho_{1,0}=0.3,\quad\rho_{2,0}=0.6,\quad\rho_{3,0}=0.1,$

while those of the commuters along the various routes are

$\bar{\rho}_{\\{1,2\\}}=0.15,\quad\bar{\rho}_{\\{1,3\\}}=0,\quad\bar{\rho}_{\\{2,3\\}}=0.29.$

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

$\bar{\rho}_{\\{1,1\\}}=0.42,\quad\bar{\rho}_{\\{2,2\\}}=0.6,\quad\bar{\rho}_{\\{3,3\\}}=0.1,$

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 (18) and the macroscopic model is (20)-(22).

Figure 5: Test 2 (commuters, Section 4.2)

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 (21),
can be computed precisely as

$\rho_{1}^{\infty}=0.255,\quad\rho_{2}^{\infty}=0.637,\quad\rho_{3}^{\infty}=0.108.$

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 (32). 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:

$\displaystyle f^{1}_{1,0}(v)=0.3\delta(v),\quad
f^{1}_{2,0}(v)=0.6\delta(v),\quad f^{1}_{3,0}(v)=0.2e^{-2v},$ $\displaystyle
f^{2}_{i,0}(v)=0,\quad\forall\,i=1,\,2,\,3.$

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

$\rho_{1,0}=\rho^{1}_{1,0}=0.3,\quad\rho_{2,0}=\rho^{1}_{2,0}=0.6,\quad\rho_{3,0}=\rho^{1}_{3,0}=0.1$

and

$\rho^{2}_{1,0}=\rho^{2}_{2,0}=\rho^{2}_{3,0}=0.$

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 (24) admits a macroscopic counterpart in closed form given by
model (30)-(31). In particular, we let

$T^{21}_{i}=0.2,\quad T^{12}_{i}=0.4,\qquad\forall\,i\in\\{1,\,2,\,3\\}.$

Figure 6: Test 3 (quarantine in the nodes, Section 4.3) with constant label
switch probabilities

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 Figures 6d, 6f) 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 [17], where the same model with quarantine is
analysed in the absence of a network. In practice, the model in [17] 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. [17, 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 [17, 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 [11] using more ad-hoc
compartmental models specifically conceived to address this phenomenon.

Figure 7: Test 3 (quarantine in the nodes, Section 4.3) with non-constant
label switch probabilities

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

$T^{21}_{i}(v)=0.2(1-e^{-v}),\quad
T^{12}_{i}(v)=0.4e^{-v},\qquad\forall\,i\in\\{1,\,2,\,3\\}.$

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., [16]. Notice that these functions are bounded from
above by the values formerly used as constant label switch probabilities:

$\sup_{v\in\mathbb{R}_{+}}T^{21}_{i}(v)=0.2,\quad\sup_{v\in\mathbb{R}_{+}}T^{12}_{i}(v)=0.4.$
(33)

For this case we have not derived a macroscopic model in closed form (though
some hints in this direction may be found in [17] in the networkless setting),
therefore in Figure 7 we only show the results of the particle model (24). 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 Figures
7d, 7f 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, Figures 7c, 7e 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 (33), 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, cf. e.g.,
[7], 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. In this
respect, our tests have essentially explored plausible prototypical scenarios,
while the systematic application of our models to real case studies will be
the object of future investigations. We remark however that 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.
[11]. 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., [16], 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.

## References

  * [1] L. Almeida, P.-A. Bliman, G. Nadin, B. Perthame, and N. Vauchelet.  Final size and convergence rate for an epidemic in heterogeneous population.  Preprint: arXiv:2010.1541, 2020.
  * [2] A. Apolloni, C. Poletto, J. J. Ramasco, P. Jensen, and V. Colizza.  Metapopulation epidemic models with heterogeneous mixing and travel behaviour.  Theor. Biol. Med. Model., 11(3):1–26, 2014.
  * [3] F. Arrigoni and A. Pugliese.  Global stability of equilibria for a metapopulation S-I-S model.  In G. Aletti, A. Micheletti, D. Morale, and M. Burger, editors, Math Everywhere, pages 229–240. Springer, 2007.
  * [4] A. D. Barbour and A. Pugliese.  Convergence of a structured metapopulation model to Levins’s model.  J. Math. Biol., 49(5):468–500, 2004.
  * [5] G. Bertaglia and L. Pareschi.  Hyperbolic models for the spread of epidemics on networks: kinetic description and numerical methods.  ESAIM Math. Model. Numer. Anal., 55(2):381–407, 2020.
  * [6] W. Boscheri, G. Dimarco, and L. Pareschi.  Modeling and simulating the spatial spread of an epidemic through multiscale kinetic transport equations.  Preprint: arXiv:2012.10101, 2020.
  * [7] M. Burger.  Network structured kinetic models of social interactions.  Preprint: arXiv:2006.15452, 2020.
  * [8] V. Colizza and A. Vespignani.  Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations.  J. Theor. Biol., 251(3):450–467, 2008.
  * [9] G. Dimarco, L. Pareschi, G. Toscani, and M. Zanella.  Wealth distribution under the spread of infectious diseases.  Phys. Rev. E, 102(2):022303, 2020.
  * [10] G. Dimarco, B. Perthame, G. Toscani, and M. Zanella.  Kinetic models for epidemic dynamics with social heterogeneity.  Preprint: arXiv:2009.01140, 2021.
  * [11] B. Espinoza, C. Castillo-Chavez, and C. Perrings.  Mobility restrictions for the control of epidemics: When do they work?  PLoS ONE, 15(7):e0235731, 2020.
  * [12] M. Garavello, K. Han, and B. Piccoli.  Models for Vehicular Traffic on Networks.  AIMS Series on Applied Mathematics. American Institute of Mathematical Sciences, Springfield, MO, 2016.
  * [13] M. Garavello and B. Piccoli.  Traffic Flow on Networks – Conservation Laws Models.  AIMS Series on Applied Mathematics. American Institute of Mathematical Sciences, Springfield, MO, 2006.
  * [14] M. J. Keeling and K. T. D. Eames.  Networks and epidemic models.  J. R. Soc. Interface, 2(4):295–307, 2005.
  * [15] W. O. Kermack and A. G. McKendrick.  Contributions to the mathematical theory of epidemics – I.  Bull. Math. Biol., 53(1–2):33–55, 1991.
  * [16] D. B. Larremore, B. Wilder, E. Lester, S. Shehata, J. M. Burke, J. A. Hay, M. Tambe, M. J. Mina, and R. Parker.  Test sensitivity is secondary to frequency and turnaround time for COVID-19 screening.  Science Advances, 7(1):eabd5393, 2021.
  * [17] N. Loy and A. Tosin.  Boltzmann-type equations for multi-agent systems with label switching.  Preprint (doi:10.13140/RG.2.2.11726.08001/2), 2021.
  * [18] M. Martcheva.  An Introduction to Mathematical Epidemiology.  Springer, 2015.
  * [19] H. Minc.  Nonnegative matrices.  Wiley-Interscience, 1988.
  * [20] L. Pareschi and G. Toscani.  Interacting Multiagent Systems: Kinetic equations and Monte Carlo methods.  Oxford University Press, 2013.
  * [21] F. Parino, L. Zino, M. Porfiri, and A. Rizzo.  Modelling and predicting the effect of social distancing and travel restrictions on COVID-19 spreading.  Preprint: arXiv:2010.05968, 2020.
  * [22] A. Wesolowski, E. zu Erbach-Schoenberg, A. J. Tatem, C. Lourenço, C. Viboud, V. Charu, N. Eagle, K. Engø-Monsen, T. Qureshi, C. O. Buckee, and C. J. E. Metcalf.  Multinational patterns of seasonal asymmetry in human movement influence infectious disease dynamics.  Nat. Commun., 8(2069):1–9, 2017.
  * [23] L. Zino and M. Cao.  Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models.  Preprint: arXiv:2013.00181, 2021.