An SIR model with viral load-dependent transmission

The viral load is known to be a chief predictor of the risk of transmission of infectious diseases. In this work, we investigate the role of the individuals’ viral load in the disease transmission by proposing a new susceptible-infectious-recovered epidemic model for the densities and mean viral loads of each compartment. To this aim, we formally derive the compartmental model from an appropriate microscopic one. Firstly, we consider a multi-agent system in which individuals are identified by the epidemiological compartment to which they belong and by their viral load. Microscopic rules describe both the switch of compartment and the evolution of the viral load. In particular, in the binary interactions between susceptible and infectious individuals, the probability for the susceptible individual to get infected depends on the viral load of the infectious individual. Then, we implement the prescribed microscopic dynamics in appropriate kinetic equations, from which the macroscopic equations for the densities and viral load momentum of the compartments are eventually derived. In the macroscopic model, the rate of disease transmission turns out to be a function of the mean viral load of the infectious population. We analytically and numerically investigate the case that the transmission rate linearly depends on the viral load, which is compared to the classical case of constant transmission rate. A qualitative analysis is performed based on stability and bifurcation theory. Finally, numerical investigations concerning the model reproduction number and the epidemic dynamics are presented.


Introduction
The route of transmission of many infectious diseases is given by social contacts among individuals.The virus shed by an infectious individual may be transmitted to a healthy one during an encounter, so that the disease also develops in the latter.There is evidence that the quantity of virus carried by the infectious individual determines the occurrence or not of the transmission: as it is reasonable to expect, higher is the viral load of the infectious individual higher is the probability of transmitting the infection.For example, the quanta emission rate (ER q ) measures the number of quanta (a quantum is the dose of airborne droplet nuclei required to cause infection in 63% of susceptible persons) into the air per time unit.The ER q for respiratory diseases (including SARS-CoV-1, SARS-CoV-2, MERS, measles virus, adenovirus, rhinovirus, coxsackievirus, seasonal influenza virus and Mycobacterium tuberculosis) has been estimated as directly [resp.inversely] proportional to the viral load in sputum [resp.the infectious dose] [23]: a more contagious strain would have higher ER q values through a higher median viral load and/or a lower infectious dose.The viral load is also the chief predictor of the risk of sexually-transmitted infections, like HIV/AIDS [25,27].
In the mathematical epidemiology community, the awareness of the importance of the viral load in the dynamics of infectious diseases has recently led to the development of epidemic models that explicitly incorporate such a microscopic trait [2,8,19,20].Specifically, in the paper [19] the authors propose a modelling framework through kinetic equations in which individuals are characterized by a discrete label and by their viral load; then, a prototype epidemic model is introduced in order to illustrate the impact of individuals' viral load on test-and-isolate activities.This work is extended in the paper [20], where the authors propose a kinetic model for the spread of an infectious disease on a graph, the nodes here representing different spatial locations.By following the wake of papers [19,20], in the paper [8] the authors introduce a compartmental susceptible-infectious-isolated-recovered model, in which the individual viral load evolves according to appropriate microscopic rules and determines the probability of isolation of infectious individuals.
To the best of our knowledge, the first epidemic model that incorporates the role of viral load in the disease transmission term has been proposed by Banerjee et al. [2].In the paper [2] an immuno-epidemiological model is introduced, where the number of susceptible people depends on the number of infectious people through the initial viral load acquired during the interactions.More precisely, according to the model [2], the growth in the number of infectious individuals increases the initial viral load, and provides a switch from the first stage of the epidemic where only people with weak immune response can be infected to the second stage where also people with strong immune response can be infected.
In the present work, we investigate how the viral load of infectious individuals affects the probability of disease transmission, and the consequent epidemic dynamics, by relying on the modelling framework of kinetic equations for multi-agent systems.Kinetic theory and Boltzmann-like equations have proved to be a very effective tool to enhance the description of infectious diseases dynamics, by allowing the incorporation in the model of not only the role of viral load [8,19,20], but also that of: social structure and wealth distribution within the host population [3,9,10,28], contact heterogeneity [11,22], implementation of lockdown measures [1] and spatial propagation of the infection [4,5,20].Specifically, we follow the approach of the papers [8,19,20], by starting from a detailed description of the microscopic dynamics of the disease spread, microscopic dynamics that is shared by all the individuals (also called the agents) that are assumed to be indistinguishable.Then, we introduce suitable kinetic equations that give a statistical portrait of the agents of the system by following exactly the prescribed microscopic rules.Eventually, from the kinetic equations we derive a macroscopic model for the aggregate description of the system that naturally inherits the details of the microscopic dynamics.
We assume that the individuals are characterized by a double microscopic state: a label, that denotes the epidemiological compartment to which they belong, and a physical quantity that is chosen to be the individual viral load.The microscopic dynamics is described in terms of microscopic interactions that allow the viral load to evolve and by the means of Markovian processes ruling the switch of compartment.We consider a basic susceptibleinfectious-recovered (SIR) compartmental structure and assume that the mechanism of disease transmission (leading the healthy individuals to become ill) depends on the viral load of the infectious individuals.
The rest of the manuscript is organized as follows.In Section 2, we present our multi-agent system and the microscopic dynamics.Then, we revise the modelling framework proposed in the paper [8] in order to derive the macroscopic compartmental model including the role of viral load in the rate of disease transmission.In Section 3, we perform a qualitative analysis of the proposed model by determining the equilibria and investigating their stability in terms of the basic reproduction number R 0 .In Section 4, some numerical simulations of the macroscopic model are performed: both the reproduction number and the epidemic temporal dynamics under the assumption of viral load-dependent rate of disease transmission are compared with those under the classical assumption of constant rate of disease transmission.Finally, in Section 5, we draw some conclusions.
2 The mathematical model: from a multi-agent system to compartmental macroscopic equations Let us consider an infectious disease spreading among individuals as a consequence of social contacts.Individuals are modelled as agents of a multi-agent system and characterized by a microscopic state.In particular, the agents are divided into disjoint compartments depending on their state of health with respect to the disease.Moreover, they are characterized by a physical quantity named viral load, that represents the quantity of viral particles present in the organism.

The microscopic model
At any time t each agent of the system is characterized by a microscopic state (x, v), where x ∈ X is a label that takes into account the epidemiological compartment to which the agent belongs, and v ∈ [0, 1] is a normalized measure of the individual's viral load, being v = 1 the maximum observable value.The evolution of both the label and the viral load may be described by the means of microscopic stochastic processes, that can be expressed through Markovian processes [8], namely through transition probabilities that is the conditional probability for an agent to change microscopic state from (j, v) to (i, v ), with (j, v), (i, v ) ∈ X × [0, 1].In general, the viral load of an individual may change both simultaneously to and independently of the switch of compartment.In the second case, we consider transition probabilities for the mere evolution of the viral load of an individual in class i ∈ X that we denote by Instead, if only the compartment changes, then we denote by P (i → j) the probability for an agent to switch from the compartment i to the compartment j.

The compartmental structure
At any time t the agents, labelled with x ∈ X , are divided in the following disjoint epidemiological compartments: • susceptible, x = S: individuals who are healthy but can contract the disease.The susceptible population increases by a net inflow, incorporating new births and immigration, and decreases due to disease transmission and natural death; • infectious, I: individuals who are infected by the disease and can transmit the virus to others.Infectious individuals arise as the result of new infections of susceptible individuals and diminish due to recovery and natural death; • recovered, x = R: individuals who have recovered from the disease after the infectious period.They come from the infectious compartments I and acquire long lasting immunity against the disease.Recovered people diminish only due to natural death.
Specifically: susceptible individuals have v ≡ 0; once infected, an individual's viral load increases until reaching a peak value (that varies from person to person) and then gradually decreases, see e.g. the representative plot of SARS-CoV-2 viral load evolution given in the paper [7], Fig. 2. Hence, for mathematical convenience [8], we assume that members of the class I are further divided into: • infectious with increasing viral load, x = I 1 ; • infectious with decreasing viral load, x = I 2 .
Note that new infections enter the class I 1 , while recovery may occur only during the stage I 2 .Finally, after the infectious period, recovered individuals may still have a positive viral load which however definitively approaches zero, as live virus could no longer be cultured (see e.g. the studies [7,17] on COVID-19 viral shedding).Also, since our model incorporates birth and death processes, we introduce the following two auxiliary compartments: individuals that enter the susceptible class by newborn or immigration, x = B, and individuals who die of natural causes, x = D. We assume that members of the class B have v ≡ 0, while those of the class D retain the viral load value at the time they died.Individuals can switch from the class B to the class S with frequency λ b and probability P (B → S) = b/ρ B (t).The quantity ρ B (t), that will be defined later, measures the size of the class B at time t.Moreover, all the living individuals can die, thus moving to the class x = D, with frequency λ µ and probability P (i → D) = µ, being i ∈ {S, I 1 , I 2 , R}.

Evolution of the viral load
Let us now focus on the mathematical modelling of the evolution of an individual's viral load v.We distinguish the two following cases when v changes over time: i) a susceptible individual, having v = 0, acquires a positive viral load (and gets infected) by interaction with an infectious individual; ii) the viral loads of infectious (I 1 , I 2 ) and recovered (R) individuals evolve naturally in virtue of physiological processes.
Given an agent labelled with S, then the necessary condition for acquiring a positive viral load is an encounter with an infectious agent (I 1 or I 2 ).Therefore, we model the disease transmission process as a binary interaction, thus relying on the typical tools of kinetic theory [24].Let us denote by λ β > 0 the frequency of these interactions.Increasing [resp.decreasing] λ β corresponds to increasing [resp.reducing] encounters among people: the lower λ β the more strengthened social distancing.
By interacting with an infectious individual carrying viral load w > 0, a susceptible individual does or does not get infected.In the first case his/her viral load after the interaction (say, v ) is positive: v > 0; in the second case it remains null: v = 0. Specifically, we consider the following microscopic binary interaction rule: where T ν β is a Bernoulli random variable of parameter ν β = ν β (w) ∈ (0, 1) describing the case of successful transmission of the disease (T ν β = 1) and the case of contact without transmission (T ν β = 0).It seems us reasonable to assume that ν β (w), that we name transmission function, is a non-decreasing function of w, the viral load of the infectious individual.
We assume that new infected individuals enter the class I 1 and they all acquire the same initial viral load, v 0 (that can be interpreted as an average initial value).We remark that this binary interaction process causes simultaneously a change of the microscopic state v and a label switch, because as soon as v becomes positive, i.e. if T ν β = 1, the susceptible individual switches to the class I 1 .In terms of transition probabilities for the susceptible individual, this can be expressed as given an encounter of the susceptible individual with an infectious one (belonging to either I 1 or I 2 ) carrying viral load w and for whom P ((i, w) Infectious and recovered individuals cannot change their viral load in binary interactions, but the evolution reflects physiological processes.Starting from the initial positive value v = v 0 , the viral load increases until reaching a given peak value and then it decreases towards zero.In this framework, the microscopic state v varies as a consequence of an autonomous process (also called interaction with a fixed background in the jargon of multi-agent systems [24]).Specifically, given an agent (I 1 , v), namely an infectious individual with increasing viral load, we consider a linear-affine expression for the microscopic rule describing the evolution of v into a new viral load v : The latter is a prototype rule describing the fact that the viral load may increase up to a certain threshold normalized to 1 by a factor proportional to (1 − v).In particular, ν 1 ∈ (0, 1) is the factor of increase of the viral load.
Similarly, given an agent (I 2 , v) or (R, v), namely an infectious individual with decreasing viral load or a recovered individual, we consider the following microscopic rule for the evolution of v: being the parameter ν 2 ∈ (0, 1) the factor of decay of the viral load.These microscopic processes happen with frequency λ γ > 0, i.e. 1/λ γ is the average increase/decay time of the viral load.
We observe here that the introduction of the sub-classes I 1 , I 2 is needed in order to implement the microscopic rules (1)-( 2) in a kinetic equation.These two rules are deliberately generic and very simple: the only aim is to distinguish individuals based on whether their viral load is increasing or decreasing and to implement two different factors ν 1 and ν 2 accordingly.
We assume that individuals in I 1 move to the class I 2 with frequency λ γ and constant probability ν 1 .In turn, individuals in I 2 move to the recovered class with frequency λ γ and constant probability ν 2 .These choices, that trace the same assumptions done in the paper [8], allow to derive λ γ ν 1 as the rate of transition from I 1 to I 2 , that is also the increase rate of the viral load.Analogously, the rate of recovery from the disease turns to be λ γ ν 2 , that is the decay rate of the viral load.Transitions I 1 → I 2 and I 2 → R are assumed to take place at the same frequency λ γ because they are driven by a common cause, namely the progression of the viral load.Hence, the rate of both transitions coincides with the progression rate of the viral load.Formally, to describe these microscopic mechanisms in terms of transition probabilities, we set

The kinetic model and the derivation of the macroscopic model
In order to give a statistical description of the multi-agent system, whose total mass is conserved in time, we introduce a distribution function for describing the statistical distribution of the individuals characterized by the pair In (3), δ(x − i) is the Dirac delta distribution centred at x = i, and f i = f i (t, v) ≥ 0 is the distribution function of the microscopic state v of the individuals that are in the ith compartment at time t.Hence, f i (t, v)dv is the proportion of individuals in the compartment i, whose microscopic state lies between v and v + dv at time t.
We assume that f (t, x, v) is a probability distribution, namely In general, the f i 's, i ∈ X , are not probability density functions because their v-integral varies in time due to the fact that individuals move from one compartment to another.
We denote by the density of individuals in the class i, thus 0 ≤ ρ i (t) ≤ 1 and Then, we define the viral load momentum of the ith compartment as the first moment of f i for each class i ∈ X , i.e.
If ρ i (t) > 0, then we can also define the mean viral load as the ratio n i (t)/ρ i (t).Instead, ρ i (t) = 0 implies necessarily f i (t, v) = 0.In this case, the mean viral load is not defined because the corresponding compartment is empty.
Starting from the microscopic dynamics illustrated in the previous section, it is possible to formally derive kinetic equations implementing exactly the microscopic processes (similarly to what done in the paper [8]).We report here the weak kinetic equations for completeness.Let ϕ : [0, 1] → R be an arbitrarily chosen test function of an observable quantity depending on the microscopic physical quantity v.For i ∈ X \ {B, D}, namely the classes of living individuals, we get: where the last term on the r.h.s.accounts for the binary interactions between susceptible individuals and infectious individuals in either , leading to the transmission of the disease, • infectious individuals with increasing viral load (i = I 1 ) • infectious individuals with decreasing viral load (i = I 2 ) • recovered individuals (i = R) As far as the disease transmission rate λ β ν β (•) is concerned, we consider that it is given by with β positive constant and p ∈ {0, 1}.
The choice p = 0 corresponds to a constant transmission function, while p = 1 reflects the experimental evidence that higher is the viral load of an infectious individual higher is his/her ability of transmitting the disease.Of course, formulations different from the linear one could be taken into account.However, since the novelty of this assumption and in absence of exhaustive field data, the linear formulation can be considered as a reasonable approximation at a first step.
In order to obtain the equations for the macroscopic densities and viral load momentum of each compartment, we set ϕ(v) = v n in ( 4)- (7), with n = 0, 1, respectively.Since we consider interaction rules that are linear in v and we assume that ν β (•) is a constant or a linear function, we obtain an exact closed system of macroscopic equations, without the need of other assumptions.This also implies that at the macroscopic level individuals in the same compartment may have heterogeneous viral loads that can be different from the mean viral load of the compartment.
The ensuing macroscopic model is given by the following system of non-linear ordinary differential equations: where we have set (with a slight abuse of notation) representing the net inflow of susceptibles and the rate of natural death, respectively.Also, for convenience of notation, in (9) we have denoted with the upper dot the time derivative and omitted the explicit dependence on time of the state variables.
From system (9) we note that the equations ruling the evolution of the densities of the compartments have an SIR structure, but the transmission term may depend on the mean viral load of the infectious population.In the simplest case that p = 0, i.e. ν β (•) is a constant function, we retrieve a classical SIR model with standard incidence [18], which reduces to ρS by noting that the differential equations for ρ R , n I1 , n I2 and n R are independent of the other ones.In such a case, the analysis of the model turns to be trivial being the equations for the densities independent of the viral load momentum.
In the present work, we take a step forward by assuming that ν β (•) is a linear increasing function, i.e. by choosing p = 1 in the system (9).With this choice, the model to be studied eventually reduces to by noting that the differential equations for ρ R and n R are independent of the other ones.
To models ( 10)- (11) we associate the following generic initial conditions Equilibria and stability properties of model (11) are investigated in the following section.
Remark 2.1.If field data concerning a specific disease showed evidence that the probability of disease transmission non-linearly depends on the viral load, one could implement a non-linear transmission function ν β (•).In such a case, it could not be possible to obtain an exact closed system of macroscopic equations, but other closure assumptions could be required.For example, in the paper [8] a monokinetic closure is used, implying that in the derivation of the macroscopic model all the individuals of a given compartment are assumed to have as viral load the mean value of that compartment.

Qualitative analysis
Let us start by ensuring that the model ( 11) is mathematically and epidemiologically well posed.It is straightforward to verify that the region with initial conditions in ( 12) is positively invariant for model (11), namely any solution of system (11) starting in D remains in D for all t ≥ 0.

The disease-free equilibrium and the basic reproduction number
The model ( 11) has a unique disease-free equilibrium (DFE), given by It is obtained by setting the r.h.s. of equations (11) to zero and considering the case ρ I1 = ρ I2 = 0.
The characteristic polynomial of J reads where with R 0 given in (14).
From the Routh-Hurwitz criterion it follows that, if R 0 < 1, then the DFE is LAS.Otherwise, if R 0 > 1, then it is unstable.
The threshold quantity R 0 is the so-called basic reproduction number for model (11), a frequently used indicator for measuring the potential spread of an infectious disease in a community.Epidemiologically, it represents the average number of secondary cases produced by one primary infection over the course of the infectious period in a fully susceptible population.
The expression of R 0 for model (11) turns out to be much more complex than that for the epidemic model ( 10) which assumes a constant disease transmission rate (we investigate more in details this point in the subsection 4.2).Note that the R 0 in ( 14) depends also on v 0 , the initial viral load of infectious individuals, a parameter that is not present in the differential equations for the densities of the compartments, namely (11a)-(11b)-(11c).
Remark 3.1.It would be interesting to investigate how the expression of the basic reproduction number R 0 varies by modifying some assumptions of model (11).For instance, one can consider the case that individuals in one between the classes I 1 and I 2 are infected but not infectious.Specifically, one can assume that • individuals in I 1 are not infectious because, for the specific disease, the period of viral load increase can be approximated to the period of latency of the infection.In such a case, the I 1 's play the role of the exposed individuals E in an SEIR model.This leads to the disappearance of the term βρ S n I1 [resp.βv 0 ρ S n I1 ] in the equation (11b) [resp.(11d)].The basic reproduction number proves to be • individuals in I 2 are not infectious because they are isolated from the community and receive treatment to decrease the viral load.This leads to the disappearance of the term βρ S n I2 [resp.βv 0 ρ S n I2 ] in the equation (11b) [resp.(11d)].In such a case, the basic reproduction number proves to be

The endemic equilibrium
Let us denote by the generic endemic equilibrium of model (11), obtained by setting the r.h.s. of equations (11) to zero and considering the case Substituting the expressions ( 16) into (11e), one gets n E I1 as a positive root of the equation Then, one can make explicit also the other components of EE: For the equilibrium to exist in D all its components must be positive.Hence, the following result can be stated.
Due to the complexity of the Jacobian matrix of system (11) evaluated at EE, we renounce to study the local stability of the endemic equilibrium.However, we make use of bifurcation analysis and show that a unique branch corresponding to the unique endemic equilibrium emerges from the criticality, namely at DFE and R 0 = 1.The emerging EE is LAS in the neighbouring of R 0 = 1 for R 0 > 1.

Central manifold analysis
To derive a sufficient condition for the occurrence of a transcritical bifurcation at R 0 = 1, we can use a bifurcation theory approach.We adopt the approach developed in [12,26], which is based on the general center manifold theory [15].In short, it establishes that the normal form representing the dynamics of the system on the central manifold is, for u sufficiently small, given by u = Au 2 + Bβu, where and Note that in (19) and ( 20) the transmission rate β has been chosen as bifurcation parameter, β is the critical value of β, x = (ρ S , ρ I1 , ρ I2 , n I1 , n I2 ) is the state variables vector, F is the r.h.s. of system (11), and z and w denote, respectively, the left and right eigenvectors corresponding to the null eigenvalue of the Jacobian matrix evaluated at criticality (i.e. at DFE and β = β).
Observe that R 0 = 1 is equivalent to so that the disease-free equilibrium is LAS if β < β, and it is unstable when β > β.
The direction of the bifurcation occurring at β = β can be derived from the sign of coefficients ( 19) and ( 20).More precisely, if A > 0 [resp.A < 0] and B > 0, then at β = β there is a backward [resp.forward] bifurcation.
For our model, we prove the following theorem.
Proof.From the proof of Theorem 1, one can verify that, when β = β (or, equivalently, when R 0 = 1), the Jacobian matrix J(DF E) admits a simple zero eigenvalue and the other eigenvalues have negative real part.Hence, the DFE is a non-hyperbolic equilibrium.
It can be easily checked that a left and a right eigenvector associated with the zero eigenvalue so that z•w = 1 are with , The coefficients A and B may be now explicitly computed.Considering only the non-zero components of the eigenvectors and computing the corresponding second derivative of F, it follows that and where z 2 , z 4 , w 4 , w 5 > 0 and w 1 < 0.Then, A < 0 < B. Namely, when β − β changes from negative to positive, the DFE changes its stability from locally asymptotically stable to unstable; correspondingly, an endemic and locally asymptotically stable equilibrium emerges.This completes the proof.

Numerical simulations
In this section, we numerically investigate how the viral load of the infectious individuals may affect the disease transmission among the population.At this aim, we compare the basic reproduction number and the numerical solutions of the macroscopic model ( 9) in the case of viral load-dependent rate of disease transmission (p = 1) with those in the classical case of constant rate of disease transmission (p = 0).Numerical simulations are performed in Matlab ® [21].We implement the 4th order Runge-Kutta method with constant step size for integrating the system (9).Platform-integrated functions are used for getting the plots.

Parametrization
Since our investigations are purely qualitative, demographic and epidemiological parameters values do not address a specific infectious disease and/or spatial area.They refer to a generic epidemic course following an SIR-like dynamics.
We are considering a model with demography and constant net inflow of susceptibles b.Since travel restrictions are usually implemented during epidemics, we assume that b accounts only for new births (which can be assumed to be approximately constant due to the short time span of our analyses).Therefore, the net inflow of susceptibles is given by where b r is the birth rate, N denotes the total resident population at the initial time, and N tot is the total (constant) system size.Note that N tot accounts for individuals belonging to all model compartments X (including B, D), whereas N refers only to living individuals.We assume an initial population of N = 10 6 individuals, representing, for example, the inhabitants of a European metropolis.The most recent data by European Statistics refer to 2020 and provide an average crude birth rate b r = 9.1/1, 000 years −1 [14] and an average crude death rate µ = 11.6/1,000 years −1 [13].The total system size N tot is set to N tot = N /(1 − bt max ), in such a way N tot = N + bt max N tot is given by the sum of the initial population, N , and the total inflow of individuals during the time interval [0, t max ], bt max N tot .The time t max is set to t max = 20 years, that is much larger than the terminal time of our numerical simulations, so ensuring that the compartment B remains not empty.
In order to make the two scenarios properly comparable, we make the following considerations.In the case S c , the quantity β c represents the rate at which infectious individuals transmit the disease in the unit of time.In the case S v , in the microscopic model the same rate is given by β v multiplied by the microscopic viral load w of the infectious individual I j , j ∈ {1, 2}; whereas, in the macroscopic model ( 9) this rate is given by β v multiplied by the mean viral load of the total infectious population: (n I1 + n I2 )/(ρ I1 + ρ I2 ).Thus, we assume that the value of β v in scenario S v is given by the value β c adopted in scenario S c rescaled by a normalization factor M ∈ (0, 1): where M represents an average quantity for (n I1 + n I2 )/(ρ I1 + ρ I2 ).It follows that For the other epidemiological parameters we take the following baseline values from the paper [8]: In particular, the product λ γ ν 1 can be interpreted as the inverse of the average time from exposure to viral load peak, whilst λ γ ν 2 as the inverse of the average time from viral load peak to recovery.All the parameters of the model as well as their baseline values are reported in Table 1.

Impact of viral load on the reproduction number
In this subsection, we investigate the impact of different modelling assumptions about the disease transmission rate on the expression and value of the reproduction number of model (9).At this aim, let us denote by the basic reproduction number of model ( 9) in the case of viral load-dependent transmission rate, S v , that is (14) with β v in place of β.It is straightforward to verify that in the case of constant disease transmission rate, S c , the basic reproduction number of model ( 9) reads (see also the paper [8]).
Remark 4.1.It is difficult to determine a priori the relationship between R v 0 and R c 0 for a given set of parameters.Nonetheless, some considerations can be made in the limit cases:  Table 1.i) λ γ 1, namely the viral load of an infected individual evolves very slowly.Then, by considering the numerator and the denominator of R v 0 and R c 0 as polynomials in λ γ and disregarding the lower order terms, one can approximate Interestingly, the ratio R v 0 /R c 0 turns to be independent of λ γ .ii) ν 1 → 0, namely all the infectious individuals have constant viral load v 0 and do not recover from the disease.
Then, R v 0 and R c 0 coincide.Indeed, iii) ν 1 , ν 2 → 1, namely in the two subsequent evolution steps after the infection, the viral load of the infected individuals reaches the maximum value 1 and then vanishes, respectively.Then, the reproduction numbers R v 0 and R c 0 read An overall view of the relationship between R v 0 and R c 0 is provided by numerically exploring their mutual position when the relevant model parameters vary in appropriate ranges.In Fig. 1, we display the relative difference of R v 0 w.r.t.R c 0 , that is the ratio as a function of the factor of viral load increase, ν 1 (Fig. 1a), the factor of viral load decay, ν 2 (Fig. 1b), and the frequency of viral load evolution, λ γ (Fig. 1c).The parameters ν 1 , ν 2 and λ γ continuously vary in the possible ranges of values  1.
We disregard the extreme cases that ν 1 , ν 2 ≈ 0 and ν 1 , ν 2 ≈ 1, which we consider to be rather unrealistic.Also, we consider nine values of the factor of transmission normalization M that span the range [0.1, 0.9], as indicated in the legend of Fig. 1.The baseline values of the other parameters are those given in Table 1.Note that the ratio R v 0 /R c 0 is independent of β c , which does not need to be assigned for the moment.From Fig. 1, we observe that the ratio R v 0 /R c 0 is a non-monotone convex function of ν 1 (Fig. 1a), an increasing function of ν 2 (Fig. 1b) and an almost constant function of λ γ (Fig. 1c), independently of M .In particular, as a function of the factor of viral load increase, ν 1 , the ratio R v 0 /R c 0 is minimum for intermediate values of ν 1 and assumes almost the same value at ν 1 = 0.05 and ν 1 = 0.95.As regards the irrelevance of the frequency of viral load evolution λ γ on the ratio R v 0 /R c 0 , it can be explained by the fact that in the current parameter setting the rate of natural death, µ, is much lower than the other parameters contributing to the reproduction numbers.From the expressions ( 22)- (23), it follows that R v 0 /R c 0 is almost independent of λ γ being the terms multiplied by µ negligible (see also the point (i) of Remark 4.1).
From Fig. 1 we also observe that the relative difference R v 0 /R c 0 − 1 decreases by increasing the factor of transmission normalization M (see (21)), eventually passing from positive to negative values (namely, from . Globally, R v 0 spans from being about 60% lower than R c 0 (M = 0.9) to 350% higher than R c 0 (M = 0.1).Also, the sensitivity of R v 0 /R c 0 to M greatly diminishes when M overcomes the value 0.5.From Fig. 1a [resp.Fig. 1b] we can note that, for a given value of M , the relative difference R v 0 /R c 0 − 1 can pass through zero by varying ν 1 [resp.ν 2 ], meaning that the mutual position between R v 0 and R c 0 changes.In particular, this happens for M = 0.4 and it is the reason why we choose it as baseline value for the next numerical investigations in this subsection.
In Fig. 2, we display the counterplots of the relative difference R v 0 /R c 0 −1 as a function of the pairs of parameters (ν 1 , ν 2 ), (ν 1 , λ γ ) and (λ γ , ν 2 ) in the ranges given in (24), by setting M = 0.4.From Fig. 2 we can evaluate the combined impact on R v 0 /R c 0 −1 of two model parameters among {ν 1 , ν 2 , λ γ } when the third one is set at the baseline value.We can note that, independently of λ γ , it is . In other words, the reproduction number of the model (9) in the scenario S v is greater than the reproduction number in the scenario S c when at least one between the factor of viral load increase and the factor of viral load decay is high, whilst R v 0 can be less than R c 0 when both ν 1 and ν 2 are medium-low.However, in any case the relative difference of R v 0 w.r.t.R c 0 is rather small: R v 0 is at most 20% greater or smaller than R c 0 .The relative difference R v 0 /R c 0 − 1 does not provide information about the exact values assumed by R v 0 and R c 0 .In order to complete the investigations, in Fig. 3 we provide the values of R v 0 and R c 0 as a function of one among the parameters {ν 1 , ν 2 , λ γ } in the ranges (24) when the other ones are set at the baseline value (Figs.3a-c), and as a function of the normalization factor M ∈ [0.1, 0.9] when the other parameters are set at the baseline value (Fig. 3d).Here, we assume that in the case of constant transmission rate the baseline value of the reproduction number is at the threshold 1, namely Other parameters values are given in Table 1.
which yields β c = 0.10 days −1 .From Fig. 3 we note that, for values of ν 1 , ν 2 and λ γ slightly smaller than the baseline value (see Figs. 3a-c, black dotted lines), it is R c 0 above the threshold 1 (blue dash-dotted lines) and R v 0 below the threshold 1 (black solid lines).This suggests that, for given epidemiological conditions, the disease dynamics predicted by model (9) can be radically different depending on the modelling assumption about the transmission function ν β (•), namely if p = 0 or p = 1 in (9).Also, from Fig. 3d, we note that the reproduction number R v 0 varies from about 3.6 to 0.4 by varying M , by crossing the threshold 1 (that is also the value of R c 0 ) for M ≈ 0.36.

Impact of viral load on the disease dynamics
In this subsection, we numerically explore the impact of different modelling assumptions about the transmission rate on the disease dynamics predicted by model (9).At this aim, we consider the following illustrative epidemiological setting.Initial data are set to the beginning of an epidemic, namely there is a single infectious individual in a totally susceptible population: Here, like in the paper [8], we assume that R c 0 = 4, which yields β c = 0.4 days −1 .
We denote by t f the terminal time of our numerical simulations (i.e.time horizon).We want that the t f is a finite time with a reasonable epidemiological interpretation.To this end, inspired by the approach adopted in the papers [6,16], we assume that t f coincides with the end of the first epidemic wave, namely t f is the first time there is less than one infectious individual in the population: In other words, t f is the first time at which ρ 1 + ρ 2 drops to 1/N tot .Of course, the presence of subsequent epidemic waves is not excluded, but for the sake of simplicity we focus here on just the first one.(26).Initial conditions and other parameters values are given in (25) and Table 1, respectively.
In order to estimate the factor of transmission normalization M , we consider the model (9) in the case S c and denote by ρ c I1 (t), ρ c I2 (t), n c I1 (t), n c I2 (t) the corresponding solutions for ρ I1 (t), ρ I2 (t), n I1 (t), n I2 (t), respectively.Then, M is set to that is the average value of the mean viral load of the infectious population over [0, t f ].In such a way, we obtain M = 0.17, yielding R v 0 = 8.47.It turns out that R v 0 is more than 110% higher than R c 0 , suggesting that the epidemic wave predicted in the scenario S v could be much more devastating than in the scenario S c .The other parameter values are given in Table 1.
In Fig. 4, we display the numerical solutions of model ( 9) in the case of viral load-dependent transmission rate, S v (black solid lines), and in the case of constant transmission rate, S c (blue dash-dotted lines).Specifically, Figs.4a and 4c [resp.4b and 4d] report the temporal dynamics of the compartment size of the infectious individuals in I 1 [resp.in I 2 ] and the corresponding mean viral load.Dotted lines indicate the terminal time (26).
In accordance with the values of the reproduction numbers, from Figs. 4a and 4b we observe that in the case S v the peak of infectious prevalence (i.e., max(ρ I1 + ρ I2 )) occurs earlier than in the case S c (at day 36 vs at day 51) and it is also higher (about 770,000 vs 600,000 total infectious individuals).However, the epidemic wave ends earlier in the case S v (t f = 176 days) than in the case S c (t f = 198 days).Cumulatively, the total number of infections during the epidemic wave are: where CI stands for cumulative incidence.We obtain that in the case S v (β = β v , p = 1) about all the initial population becomes infected, whilst in the case S c (β = β c , p = 0) only 6,000 infections are avoided.From a mathematical point of view, this means that the area under the curve ρ I1 + ρ I2 changes little by varying the modelling assumption about the disease transmission rate.In In Figs.4c and 4d we report the temporal dynamics of the mean viral load of the infectious compartments.As anticipated in Section 2, we highlight that the mean viral load is reliable when the number of particles is sufficiently high due to the law of large numbers (see [8] for a more detailed discussion).From Fig. 4c we note that the mean viral load of the compartment I 1 in the case S v (black solid line) is -for most of the time horizon -higher than in the case S c (blue dash-dotted line).At variance, the mean viral load of the compartment I 2 is smaller in the case S v w.r.t. the case S c (Fig. 4d).This suggests that the model ( 9) under the assumption of constant disease transmission tends to underestimate the mean viral load of infectious individuals in the increasing phase and to overestimate the mean viral load of infectious individuals in the decreasing phase.From Figs. 4c and 4d we also note that, at the end of the time horizon when the compartments I 1 and I 2 are almost empty, the mean viral load remains approximately constant at a positive value, suggesting that the viral load momentum n I1 [resp.n I2 ] and the density ρ I1 [resp.ρ I2 ] go to zero with the same speed.

Conclusion
In this work, we propose and analyse an SIR epidemic model with viral load-dependent transmission.The compartmental model is formally derived -by the means of kinetic equations -from a stochastic particle description of the individual course of the disease and the viral load progression.This approach allows the macroscopic model to inherit the features of the microscopic dynamics related to the heterogeneity of the viral load in the population [8].
The main results are as follows: • the particle stochastic model provides that, in the binary interaction between a susceptible and an infectious individual, the probability for the former to get infected depends on the viral load of the latter.In particular, the transmission function is a non-decreasing function of the viral load of the infectious individual.In the macroscopic model, the rate of disease transmission turns out to be a function of the mean viral load of the infectious population; • we analytically and numerically investigate the impact of different modelling assumptions about the disease transmission rate on the epidemic dynamics.In particular, we consider the case that the transmission rate linearly depends on the viral load (scenario S v ), which is compared to the classical case of constant transmission rate (scenario S c ); • we determine explicitly the equilibria of the macroscopic model in the scenario S v and study their stability in terms of the basic reproduction number (R 0 ) of the model.We prove that a transcritical forward bifurcation occurs at the disease-free equilibrium and R 0 = 1.The expression of the R 0 appears to be much more complex than that in the scenario S c , by depending -inter alia -on the initial viral load of the infectious individuals; • the numerical simulations unravel the relationship between the reproduction numbers R 0 in the scenarios S v and S c when the model parameters vary in appropriate ranges.We observe that the mutual position between the two R 0 's is almost independent of the frequency of the viral load evolution, but it may be noticeably affected by the factors of viral load increase and decay.Interestingly, for given parameter values, it may be R 0 < 1 in scenario S v and R 0 > 1 in scenario S c , suggesting that the disease dynamics can be radically different depending on the modelling assumption about the transmission rate; • we simulate an epidemic wave by assuming R 0 = 4 in the scenario S c and estimating the R 0 in the scenario S v accordingly.We obtain that in the case of viral load-dependent transmission, the epidemic wave is more severe, with a higher and earlier prevalence peak, than in the case of constant transmission.Also, the model in the scenario S c tends to underestimate the mean viral load of infectious individuals whose viral load is increasing and to overestimate the mean viral load of infectious individuals whose viral load is decreasing.
The role of viral load in the dynamics of infectious diseases has recently attracted the interest of mathematical epidemiologists.Our work makes a further step in this line of research, by serving as a proof-of-principle verification of the impact of the individuals' viral load on the disease transmission rate and the consequent epidemic dynamics.
In the proposed framework, the description of the microscopic mechanisms and the heterogeneity of the viral load at the microscopic level allows one to derive a macroscopic model, which provides for a richer description of the disease spreading in the host population w.r.t.classical epidemic models.Here we only consider the explicit influence of the viral load on the transmission mechanism, but, in principle, other switches of individuals between compartments may depend on the viral load at the microscopic level, and on the mean viral load at the macroscopic level.Also, more complex situations could be addressed, for example by assuming different initial viral loads of the infectious individuals that may give rise to a different epidemic scenario.
We underline that our model does not address a specific infectious disease.Of course, in presence of exhaustive data concerning the progression of the individuals' viral load, the modelling assumptions can be reformulated or adjusted according on the particular case.However, even if our model is too simple to provide reliable solutions to real-world epidemics from the quantitative point of view, our theoretical findings can be used to inform more complex simulation models developed for specific epidemiological scenarios where more realistic descriptions of the biological and epidemiological processes are included.

Figure 1 :
Figure 1: Relative difference of the reproduction number of the model (9) in the scenario S v , R v 0 , w.r.t. the reproduction number in the scenario S c , R c 0 , for nine values of the factor of transmission normalization, M ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9} (as indicated in the legend).Panel (a): R v 0 /R c 0 − 1 as a function of the factor of viral load increase, ν 1 .Panel (b): R v 0 /R c 0 − 1 as a function of the factor of viral load decay, ν 2 .Panel (c): R v 0 /R c 0 − 1 as a function of the frequency of viral load evolution, λ γ .Black dotted lines indicate the corresponding baseline value.Other parameters values are given in Table1.

Figure 2 :
Figure 2: Counterplots of the relative difference of the reproduction number of the model (9) in the scenario S v , R v 0 , w.r.t. the reproduction number in the scenario S c , R c 0 .Panel (a): R v 0 /R c 0 − 1 versus the factor of viral load increase, ν 1 , and the factor of viral load decay, ν 2 .Panel (b): R v 0 /R c 0 − 1 versus the factor of viral load increase, ν 1 , and the frequency of viral load evolution, λ γ .Panel (c): R v 0 /R c 0 − 1 versus the frequency of viral load evolution, λ γ , and the factor of viral load decay, ν 2 .The intersection between white dotted lines indicates the corresponding baseline value.Other parameters values are given in Table1.

Figure 3 :
Figure 3: Reproduction number of the model (9) in the scenario S v , R v 0 (black solid lines), and reproduction number in the scenario S c , R c 0 (blue dash-dotted lines).Panel (a): R v 0 and R c 0 as functions of the factor of viral load increase, ν 1 .Panel (b): R v 0 and R c 0 as functions of the factor of viral load decay, ν 2 .Panel (c): R v 0 and R c 0 as functions of the frequency of viral load evolution, λ γ .Panel (d): R v 0 and R c 0 as functions of the factor of transmission normalization, M .Black dotted lines indicate the corresponding baseline values.Red x-marks indicate the intersection points between R v 0 and R c 0 .Other parameters values are given in Table1.

Figure 4 :
Figure 4: Numerical solutions as predicted by the model (9) in scenarios S v (black solid lines) and S c (blue dashdotted lines).Panel (a): compartment size of infectious individuals with increasing viral load, I 1 .Panel (b): compartment size of infectious individuals with decreasing viral load, I 2 .Panel (c): mean viral load of infectious individuals with increasing viral load, I 1 .Panel (d): mean viral load of infectious individuals with decreasing viral load, I 2 .Dotted lines indicate the corresponding terminal time, as defined in(26).Initial conditions and other parameters values are given in(25) and Table1, respectively.

Table 1 :
List of model parameters with corresponding description and baseline value.

Table 2 :
(25)vant quantities as predicted by the model(9)in the case of viral load-dependent transmission rate (scenario S v , first line) and in the case of constant transmission rate (scenario S c , second line).First column: infectious prevalence peak, max(ρ I1 +ρ I2 )N tot .Second column: time of infectious prevalence peak, argmax(ρ I1 +ρ I2 ).Third column: terminal time, t f .Fourth column: cumulative incidence at t f , CI.Initial conditions and other parameters values are given in(25)and Table1, respectively.