key: cord-0906805-ksj6racy
authors: Loy, Nadia; Tosin, Andrea
title: A viral load-based model for epidemic spread on spatial networks
date: 2021-04-25
journal: Mathematical biosciences and engineering : MBE
DOI: 10.3934/mbe.2021285
sha: 4d7d6aec60b20538d41ecf1db9226621f9f3121e
doc_id: 906805
cord_uid: ksj6racy

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.

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 Boltzmanntype 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 Boltzmanntype 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 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 ∈ I, where I is a finite ordered index set with |I| = n ∈ N, for instance I = {1, . . . , n} ⊂ 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 P := [P ij ] i,j∈I ∈ R n×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 P is irreducible.

Besides their current vertex, individuals are further characterised by a microscopic state v ∈ 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 = pv + qv * , v * = pv * + qv,

where v, v * represent the pre-interaction viral loads of two interacting individuals and v , v * their post-interaction viral loads. Moreover, p, q ∈ 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 * , 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 ∈ I, from node to node. In particular, we fix

where ν 1 ∈ [0, 1] is the physiological decay rate of the viral load of an infected individual, η is a random coefficient modelling stochastic fluctuations of the physiological decay rate and ν 2 ∈ [0, 1] is the transmission rate of the contagion. Denoting by · the expectation with respect to the law of η, we assume η = 0 and moreover η ≥ ν 1 − 1 so as to meet the condition p ∈ R + , which, together with q ∈ R + , guarantees v , v * ≥ 0 for all v, v * ≥ 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 ∈ I in the network and its viral load V t ∈ R + at time t. Then we describe the evolution of X t , V t during a time step ∆t > 0 by means of the following update rules:

where:

i) Θ, Ξ ∈ {0, 1} are independent Bernoulli random variables discriminating whether a vertex jump or a change in the viral load take place (Θ, Ξ = 1) or not (Θ, Ξ = 0) in the time step ∆t. In particular, we let Prob(Θ = 1) = χ∆t, Prob(Ξ = 1) = µ∆t,

χ, µ > 0 being the mobility and social contact rates, and we assume ∆t ≤ min{ 1 χ , 1 µ } for consistency;

ii) J t ∈ I is a random variable returning the new vertex after a vertex jump, with

iii) V t ∈ R + is the new viral load after a social contact. In view of (3) we have V t = pV t + qV * t , where V * t denotes the viral load at time t of the other individual participating in the interaction; iv) δ Xt,X * 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 * t ∈ 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 ∈ R + at time t > 0 in vertex i ∈ I, that we denote by f i = f i (v, t) : R + × R + → R + . Taking advantage of the procedure detailed in [17] , see also [20] for further reference, in the continuoustime limit ∆t → 0 + we find that f i satisfies formally the following Boltzmann-type kinetic equation in weak form:

where ϕ : R + → 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 (3)-(4) as a function of the integration variables v, v * .

We observe that (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 ϕ(v) = 1 in (7) we obtain that ρ i satisfies the equation

which in vector notation reads dρ dt = χ(P − I)ρ

with ρ := (ρ i ) i∈I . From (10) we can investigate the stationary mass distribution ρ ∞ ∈ R n + emerging for large times, namely the vector satisfying the equation

In practice, ρ ∞ should be an eigenvector of P corresponding to the eigenvalue 1. Notice that P admits indeed this eigenvalue, which coincides also with its spectral radius, because it is a stochastic matrix. Moreover, since 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 ρ ∞ , but we need to manage its non-uniqueness (scalar multiples are also eigenvectors). For this, we notice that in our case the 1 -norm of ρ, i.e. the quantity ρ 1 := i∈I ρ i , is conserved in time, indeed from (9) and (2) we deduce d dt

Thus ρ ∞ 1 is fixed by the initial condition, which removes the remaining degree of freedom in the identification of a unique physically admissible ρ ∞ ∈ R n + solving (11) .

Still from Perron-Frobenius theory we have that 1 is also the spectral limit of P, namely the maximum real part of its eigenvalues. Therefore the real parts of the eigenvalues of the matrix P − I are non-positive and, in particular, the unique eigenvalue with null real part is simple. This implies that ρ ∞ 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 1 -norm discussed above.

Summarising, we have proved:

There exists a unique physically admissible solution ρ ∞ ∈ R n + to (11) , which is a stable and attractive asymptotic density distribution for (10) . (7) we may investigate the evolution of the mean viral load

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

From here, invoking (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 m ∞ = 0, which represents the eradication of the infection in all nodes of the network. Linearising the previous equation around the equilibrium

Next, recalling that ρ ∞ has strictly positive components we set c := min i∈I ρ ∞ i > 0 and we use

which for t → +∞ implies i∈I ρ ∞ i m i → 0 and therefore m i → 0 for all i ∈ I. In this case the configuration m ∞ = 0 is locally asymptotically stable, hence the infection may be globally eradicated in the long run;

which for t → +∞ yields i∈I ρ ∞ i m i → +∞ and therefore m i → +∞ for some i ∈ I. In this case the configuration m ∞ = 0 is unstable: the infection is not under control and tends to grow in time;

iii) if ν 1 = ν 2 then i∈I ρ i m i is constant in time. This indicates that the configuration m ∞ = 0 is stable but not attractive: the infection is under control but cannot be eradicated and becomes endemic.

If the decay and transmission rates ν 1 , ν 2 are not constant but vary from vertex to vertex of the graph, i.e. ν 1 = ν 1,i and ν 2 = ν 2,i , then (13) modifies as

Apart from the relatively trivial cases in which either ν 1,i > ν 2,i , ν 1,i < ν 2,i or ν 1,i = ν 2,i for all i ∈ 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

where 0 < 1 is a small parameter playing conceptually the role of the Knudsen number of classical kinetic theory. In the hydrodynamic limit → 0 + a splitting of (7) is possible in: i) intra-vertex interactions without jumps on the quick time scale τ := t/ , ruled by

which lead f i to converge for large τ to a local equilibrium distribution M i (the local Maxwellian) parametrised by the macroscopic quantities conserved by the interactions (the collisional invariants) in the ith 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 intravertex interactions.

Taking ϕ(v) = 1 in (14) we discover that in all vertices the density ρ i is constant on the time scale τ . Taking instead ϕ(v) = v we get

Let us assume first that

Plugging these local Maxwellians into (15) and choosing ϕ(v) = 1 yields the same time evolution of the densities ρ i , i ∈ I, as (9) . Then, in particular, ρ i → ρ ∞ i > 0 for all i ∈ I when t → +∞ as asserted by Proposition 2.1. Choosing instead ϕ(v) = v and i = i * we may investigate the evolution of the mean viral load in the i * 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 * i * < 1, whence ρ i * m i * → 0 for t → +∞. Considering that ρ i * converges to a non-zero asymptotic value, we finally infer m i * → 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 τ -scale. Next, let us consider the opposite situation, namely ν 1,i = ν 2,i for all i ∈ I \ {i * } and ν 1,i * > ν 2,i * . In this case, in all vertices i = i * the local contagion dynamics yield local Maxwellians M ρi,mi (v) parametrised by both ρ i and m i while in the i * th vertex it results m i * → 0 as τ → +∞, thus the local Maxwellian is

where we have used the fact that i∈I\{i * } P ij = 1 − P i * j . Let us introduce the set of indices of the vertices different from vertex i * and directly linked to the latter:

Notice that I * is non-empty due to the strong connectivity of the graph (if P i * j = 0 for all j = i * then vertex i * could not be reached from any other vertex). Let a := min j∈I * P i * j > 0, then from (17) we deduce d dt

whence, integrating both sides in time and considering that i∈I\{i

By Grönwall's inequality we deduce then

which says that i∈I * ρ i m i → 0, hence that m i → 0 for all i ∈ I * as t → +∞. 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 I * = I \ {i * }, i.e. if all vertices of the graph are directly linked to vertex i * , 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 ν 1 = ν 2 . Finally, if ν 1,i * < ν 2,i * for some i * ∈ I then from (16) we easily compute that the intra-vertex interactions produce locally a blow-up of the mean viral load: m i * → +∞ as τ → +∞. The mean viral load will consequently blow in every other vertex i ∈ I where ν 1,i ≤ ν 2,i while the infection will be eradicated only in the vertices where ν 1,i > ν 2,i .

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 ) ∈ I × 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 ∈ R + be their viral load. We model the evolution of (X o,t , X d,t ) and V t during a time step ∆t > 0 as

where Θ, Ξ are like in (6), V t = pV t + qV * t from (3) and V * t , X * 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 ) ∈ I × 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 ,Xo,t ∈ [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 (18) an aggregate description in terms of the distribution function of the viral load v ∈ R + of commuters between vertices i, j ∈ 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) : R + × R + → R + such a distribution function. Using again the technique illustrated in [17] with minimal adaptations, in the continuous-time limit ∆t → 0 + we formally obtain from (18) that f ji satisfies the following Boltzmann-type kinetic equation in weak form:

where v is given by (3)- (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 (19) features two important differences with respect to (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 ijcommuters 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 = 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

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

Notice that the time evolution of ρ ji depends only on ρ ji itself and on ρ 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 ρ ij we easily see that the total density of ij-commuters, namely ρ ji + ρ ij , is constant in time. Let us denote by 1ρ {i,j} := ρ ji + ρ ij the total mass of ij-commuters. Then ρ ij =ρ {i,j} − ρ ji , whence

which implies that ρ ji , ρ 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 (20) , the evolution equation for the total density ρ i turns out to be

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 ∈ I and this equation becomes

which reminds more closely of (9) but with P ij ρ j replaced by ρ ij . Notice that, unlike (9), this is not a self-consistent equation in terms of the ρ 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 (19) with ϕ(v) = v and taking also (3)- (4) and (20) 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 (22)). Contagion may spread along all routes because of mixed social contacts in each vertex.

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 ∈ 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 → 2 means that the individual is quarantined; conversely, the label switching 2 → 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 ∈ I to which they belong, the label H t ∈ {1, 2} identifying their quarantine state and their viral load V t ∈ R + . Next, we model the state update in a time step ∆t > 0 as

where Λ ∈ {0, 1} is a Bernoulli random variable such that Prob(Λ = 1) = λ∆t which discriminates whether a label switching takes place (Λ = 1) or not (Λ = 0) in the time step ∆t and λ > 0 is the rate of viral testing. Moreover, I t ∈ {1, 2} is the new label assigned to the individual after a viral test, with in particular

The random variables Θ, Ξ are like in Section 2.1 but in this case we fix ∆t ≤ min{ 1 χ , 1 λ , 1 µ } 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 (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 (3)-(4) into a rule expressing a progressive reduction of viral load due to the lack of social contacts:

which is of the form (3) with p = 1 − γ + ξ and q = 0. Here, γ ∈ [0, 1] is the decay rate of the viral load of a quarantined individual, which in general is expected to satisfy γ ≥ ν 1 because quarantined individuals may be exposed to specific medical treatments in addition to the physiological recovery. Moreover, ξ ∈ R is a centred random coefficient modelling stochastic fluctuations of the decay rate.

With the restriction ξ ≥ γ − 1 we guarantee v ≥ 0 for all v ≥ 0. On the whole, we set

where H * t ∈ {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) : R + × R + → R + be the distribution function of the viral load v ∈ R + at time t > 0 in the vertex i ∈ I for an individual with label h ∈ {1, 2}. Combining (7) and the results in [17] we deduce from (24) that in the continuous-time limit ∆t → 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 (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 * , t). Explicitly, the equations for the two compartments h = 1, 2 read then:

• non-quarantined individuals (h = 1):

with v given by (3)-(4);

• quarantined individuals (h = 2):

with v given by (25).

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 

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

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 ρ i 's alone. This is the case of the model with quarantine, as we have in general P 1 ij = 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 ρ h i :

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

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

Letting ϕ(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),

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 |I| = 1) for both constant and non-constant label switch probabilities.

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. Figure 3 ) ( Figure 5) (Figures 6, 7 ) In all tests we consider the strongly connected graph illustrated in Figure 2 , which models a network of e.g., three different cities. Hence I = {1, 2, 3}. We choose the transition matrix as P =   P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 Table 1 .

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

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: 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 ν 1 < ν 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 ≥ 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 ν 1 < ν 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 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 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.

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 (32). 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 → 2 and 3 → 3, where m 32,0 = m 33,0 = 1. Therefore node 3 is again the infection hotspot. As a matter of fact, 3 → 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 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

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 |ρ ∞ i − ρ i,0 | 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.

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 nonquarantined individuals coinciding with the corresponding P ij 's in (32). The transition probabilities of quarantined individuals are instead the trivial ones P 2 ij = δ 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 ρ 1,0 = ρ 1 1,0 = 0.3, ρ 2,0 = ρ 1 2,0 = 0.6, ρ 3,0 = ρ 1 3,0 = 0.1 and ρ 2 1,0 = ρ 2 2,0 = ρ 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 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 → +∞. 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.

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., [16] . 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 [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 

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.

Final size and convergence rate for an epidemic in heterogeneous population

Metapopulation epidemic models with heterogeneous mixing and travel behaviour

Global stability of equilibria for a metapopulation S-I-S model

Convergence of a structured metapopulation model to Levins's model

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

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

Network structured kinetic models of social interactions

Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations

Wealth distribution under the spread of infectious diseases

Kinetic models for epidemic dynamics with social heterogeneity

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

Models for Vehicular Traffic on Networks

Traffic Flow on Networks -Conservation Laws Models

Networks and epidemic models

Contributions to the mathematical theory of epidemics -I

Test sensitivity is secondary to frequency and turnaround time for COVID-19 screening

Boltzmann-type equations for multi-agent systems with label switching

An Introduction to Mathematical Epidemiology

Nonnegative matrices

Interacting Multiagent Systems: Kinetic equations and Monte Carlo methods

Modelling and predicting the effect of social distancing and travel restrictions on COVID-19 spreading

Multinational patterns of seasonal asymmetry in human movement influence infectious disease dynamics

Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models

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 evolution-ary 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.