Cytomegalovirus dynamics model with random behavior

1.   Introduction

Mathematical models of infectious diseases spreading have played a significant role in infection control. On one hand, they have given an important contribution to the biological and epidemiological understanding of disease outbreak patterns; on the other hand, they have helped to determine how and when to apply control measures in order to quickly and most effectively contain epidemics [1]. Research in this field is constantly evolving and ever new challenges are launched from the real world (just think of the ongoing COVID–19 pandemic). One among the many increasingly attractive topics is the mutual influence between the individual behaviours and choices and the disease dynamics [2,30].

In mathematical epidemiology literature a prominent position is occupied by the compartmental epidemic models. They are macroscopic models where the total population is divided into disjoint compartments according to the individual status with respect to the disease, and the switches from a compartment to another follow given transition rules. The size of each compartment represents a state variable of the model, whose rate of change is ruled by a balance differential equation. The milestone of compartmental models is the well–known deterministic Susceptible–Infected–Removed (SIR) model, proposed by Kermack and McKendrick in 1927 [18].

Like any mathematical model, also epidemic models postulate some simplifying assumptions that are needed to make them analytically tractable and/or numerically solvable. Quantifying the impact of such simplifications is extremely important to understand the model reliability and identify its range of application. For example, deterministic compartmental models are valid for large populations. Hence, they can hardly describe situations in which compartments are almost empty (for example at the onset of an epidemic, when the infectious individuals are very few) and, then, stochastic fluctuations cannot be disregarded.

A significant aspect neglected by classical epidemic models is the heterogeneity of disease transmission and progression linked to the viral load of each infected individual. Viral load is defined as a quantitative viral titre (e.g. copy number) [11] and may represent a useful marker for assessing viral kinetics, disease severity and prognosis. Indeed, symptoms and mortality induced by the infection may depend on the individual viral load, like asserted, for example, by studies on seasonal flu [21], measles [28] and COVID–19 disease [14]. The quantity of virus in the organism can also influence the results by screening and diagnostic tests, which are capable of detecting a different quantity of virus per $ml$ according to their sensitivity. Hence, the viral load affects the probability for an individual of being diagnosed and, consequently, home isolated or hospitalized, thus preventing the possibility of him/her infecting other people. In this context, assessing the interplay between the frequency of testing and sensitivity of the tests is crucial for planning prevention and mitigation measures [20]. Also the timing of testing is fundamental: for example in the case of acute rubella, in order to have laboratory confirmation of infection, viral specimen should be collected as soon after symptom onset as possible, preferably one to three days after onset, but no later than seven days post–onset [5]. Last but not least, the viral load can be a strong determinant of transmission risk [15], and the knowledge of the duration of viral shedding plays a key role in tracing the evolution of the infectious disease [6]. For example, it is estimated that SARS–CoV–2 viral load peaks just before the symptom onset, i.e. during the pre–symptomatic stage of infection [17,11], and pre–symptomatic patients are responsible of about $44\%$ of secondary infections [17]. In the case of congenital rubella syndrome, infants can shed the virus up to one year, but samples should be collected prior to three months of age because by three months of age approximately 50% will no longer shed virus [5].

The mathematical framework of multi–agent systems [27] allows one to introduce a detailed microscopic description of the interactions between individuals, that are generally called agents, within a population. One of the key aspects is that it allows one to recover a statistical description of the system by introducing a probability density function accounting for the statistical distribution of some microscopic traits of the individuals. Its evolution may be described by kinetic equations that also permit to derive macroscopic equations, i.e. macroscopic models, that, thus, inherit a large number of features of the original microscopic dynamics. In particular, the authors in [22] introduced a particle model describing a microscopic dynamics in which agents have a double microscopic state: a discrete label that switches as a consequence of a Markovian process and a microscopic trait that changes as a consequence of binary interactions or interactions with a background. The trait may take into account the individual viral load, while the label denotes the compartment to which the agent belongs. The authors then derived nonconservative kinetic equations describing the evolution of the distribution of the microscopic trait for each label and, eventually, macroscopic equations for the densities and momentum of the microscopic trait of each compartment.

Kinetic equations have been applied to compartmental epidemic models in order to take into account the role of wealth distribution during the spread of infectious diseases, for example in [8,9]. In these works, the authors described in more detail social contacts among the individuals but still relied on an SIR–like structure to model contagion dynamics. To our best knowledge, the only kinetic model taking into account the spread of an infectious disease by means of interactions and including the individuals' viral load is the one proposed in [25], where, however, the authors did not consider epidemiological compartments.

Motivated by the previous arguments, in the present work we propose a microscopic stochastic model allowing one to describe the spread of an infectious disease as a consequence of the interactions among individuals who are characterized by means of their viral load. Once infected, the viral load of the individuals increases up to a maximum peak and then decreases as a consequence of a physiological process. Furthermore, the individuals are labelled in order to indicate their belonging to one of the disjoint epidemiological compartments. Specifically, we consider an SIR–like dynamics with an isolation mechanism that depends on the individual viral load (Section 2). We derive kinetic evolution equations for the distribution functions of the viral load of the individuals in each compartment and, eventually, a macroscopic model for the densities and viral load momentum (Section 3). We perform then a qualitative analysis of the ensuing macroscopic model (Section 4), and we present some numerical tests of both the microscopic and the macroscopic models to show the matching between the aggregate trends obtained from the macroscopic descriptions and the original particle dynamics simulated by a Monte Carlo approach (Section 5). Finally, we draw some conclusions and we briefly sketch possible research developments (Section 6).

2.   A multi–agent system describing the disease spread through interactions

Let us consider a large system of interacting individuals in presence of an infectious disease that spreads through social contacts. The total population at time $t$ is divided into disjoint epidemiological compartments according to the health status with respect to the disease, to each of which we associate a label $x \in \mathcal{X}$. Individuals, that we shall also call the agents, are characterized by the evolution stage of the disease–related viral load, that is the number of viral particles present in the organism. Let us denote with $v\in[0, 1]$ a normalized measure of the individual viral load at time $t$, where $v = 1$ represents the maximum observable value. We want to describe the microscopic mechanisms modelling the interactions between individuals, which are the means of the transmission of the disease, and the switch of compartment of each individual that follows from the disease progression.

Being our final aim the proposal of a macroscopic model, we derive as an intermediate stage a statistical description of our multi–agent system through kinetic equations, by which we then derive macroscopic equations. 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 agents characterized by the pair $(x, v) \in \mathcal{X}\times[0, 1]$, as

where $\delta(x-i)$ is the Dirac delta distribution centred at $x = i$, and we assume that $f(t, x, v)$ is a probability distribution, namely

In (1)–(2), $f_i = f_i(t, v)\geq 0$ is the distribution function of the microscopic state $v$ of the agents that are in the $i$th compartment at time $t$. Hence, $f_i(t, v)dv$ is the proportion of agents in the compartment $i$, whose microscopic state lies between $v$ and $v+dv$ at time $t$. In general, the $f_i$'s, $i \in \mathcal{X}$, are not probability density functions because their $v$–integral varies in time due to the fact that agents move from one compartment to another. We denote by

the density of agents in the class $i$, thus $0\leq\rho_i(t)\leq 1$ and

Then, we define the viral load momentum of the $i$th compartment as the first moment of $f_i$ for each class $i \in \mathcal{X}$, i.e.

If $\rho_i(t)>0$, then we can also define the mean viral load as the ratio ${n_i(t)}/{\rho_i(t)}.$ We observe that $\rho_i(t) = 0$ implies instead necessarily $f_i(t, v) = 0$ and, therefore, $n_i(t) = 0$. In such a case, the mean viral load is not defined because the corresponding compartment is empty. We also remark that, if the compartment is almost empty, then the mean viral load $n_i/\rho_i$, $i\in \mathcal{X}$, might not be fully consistent with the empirical mean viral load resulting from the particle description because the law of large numbers does not apply.

2.1.   The compartmental structure

The individuals, labelled with $x\in \mathcal{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. We assume that members of this class are asymptomatic or mildly symptomatic, hence they move freely. Infectious individuals arise as the result of new infections of susceptible individuals and diminish due to recovery and natural death or because they are identified and isolated from the general population;

● isolated, $H$: infected individuals who have been identified and fully isolated from the general population by home isolation or hospitalization. Members of this class come from the infectious compartment $I$ and get out due to recovery or death. We assume that this class includes patients showing severe symptoms, that can also die due to the disease;

● recovered, $x = R$: individuals who have recovered from the disease after the infectious period. They come from the infected compartments $I$ and $H$ and acquire long lasting immunity against the disease.

Specifically: susceptible individuals have $v\equiv0$; 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 [6], Fig. 2. Hence, for mathematical convenience, we assume that members of classes $I$ and $H$ are further divided into:

● infectious, $x = I_1$, and isolated, $x = H_1$, with increasing viral load;

● infectious, $x = I_2$, and isolated, $x = H_2$, with decreasing viral load.

Note that new infections enter the class $I_1$, while recovery may occur only during the stages $I_2$ or $H_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 [6,17] on COVID–19 viral shedding).

Also, since our model incorporates birth and death processes, we introduce the following three auxiliary compartments: individuals that enter the susceptible class by newborn or immigration, $x = B$; individuals who die of natural causes, $x = D_\mu$; and individuals who die from the disease, $x = D_d$. We assume that members of class $B$ have $v\equiv 0$, while those of classes $D_\mu$ and $D_d$ retain the viral load value at time they died.

2.2.   Evolution of the viral load

Let us now focus on the mathematical modelling of the evolution of an individual viral load $v$. We distinguish the two following cases when $v$ changes over time: ⅰ) a susceptible individual, having $v = 0$, acquires a positive viral load (and get infected) by interaction with an infectious individual; ⅱ) the viral loads of infected ($I_1$, $I_2$, $H_1$, $H_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 individual ($I_1$ or $I_2$). Let us denote with $\lambda_\beta>0$ the frequency of these interactions. Increasing [resp. decreasing] $\lambda_\beta$ corresponds to increasing [resp. reducing] encounters among people: the lower $\lambda_\beta$ the more strengthened social distancing.

By interacting with an infectious individual, a susceptible individual may or may 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 rule:

where $T_{\nu_\beta}$ is a Bernoulli random variable of parameter $\nu_\beta\in (0, 1)$ describing the case of successful contagion when $T_{\nu_\beta} = 1$ and the case of contact without contagion when $T_{\nu_\beta} = 0$. 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_{\nu_\beta} = 1$, the susceptible individual switches to the class $I_1$.

Infectious, isolated and recovered individuals cannot change their viral load in binary interactions, but the evolution reflects physiological processes. In particular, 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. The peak can be reached in either the stage $I_1$ or $H_1$, i.e. before or after an infectious individual is possibly isolated.

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 [27]). Specifically, given an agent $(I_1, v)$ or $(H_1, v)$, namely an infected 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, $\nu_1\in (0, 1)$ is the factor of increase of the viral load.

Similarly, given an agent $(I_2, v)$, $(H_2, v)$ or $(R, v)$, namely an infected individual with decreasing viral load or a recovered individual, we consider the following microscopic rule for the evolution of $v$:

being the parameter $\nu_2\in(0, 1)$ the factor of decay of the viral load. These microscopic processes happen with frequency $\lambda_\gamma>0$. We observe here that the introduction of the sub–classes $I_1, \, I_2$ and $H_1, \, H_2$ is needed in order to implement the microscopic rules (3)–(4) 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 $\nu_1$ and $\nu_2$ accordingly.

2.3.   A microscopic stochastic model

Let us now define a microscopic stochastic process implementing the modelling assumptions defined so far. Let us consider an agent characterized by the pair of random variables $(X_t, V_t)$, where $X_t \in \mathcal{X}$ is the label denoting the compartment to which the agent belongs and $V_t\in [0, 1]$ is the viral load. Then, the random variable $X_t$ changes in consequence of a Markovian jump process, while the microscopic state $V_t$ may change either because of a binary interaction with an agent characterized by $(Y_t, W_t)$, $Y_t \in \mathcal{X}$, $W_t \in [0, 1]$ or because of an autonomous process according to the discussion in Section 2.2. In particular, two different types of stochastic processes may happen:

$1.$ the agents may change their label according to a process that is independent of the change of the viral load: birth and death processes, isolation process for individuals whose viral load is increasing $(I_1\rightarrow H_1)$, isolation process for individuals whose viral load is decreasing $(I_2\rightarrow H_2)$, where the notation $(j\rightarrow i)$ indicates the switch from compartment $j$ to compartment $i$;

$2.$ the agents may change both their viral load and their label simultaneously: $(S\rightarrow I_1)$, $(I_1\rightarrow I_2)$, $(I_2\rightarrow R)$, $(H_1\rightarrow H_2)$, $(H_2\rightarrow R)$. This class of processes also includes the evolution of the viral load of individuals who remain in the same compartment, i.e. $(i\rightarrow i)$, $\forall i \in \mathcal{X}$.

The two stochastic processes above may be expressed in the following rule describing the variation of $X_t$ and $V_t$ of a generic representative agent of the system during a time interval $\Delta t>0$:

where $\Sigma$ and $\Psi$ are indicator functions. In particular, $\Sigma = 1$ if the agents of the compartment labelled with $X_t$ change both their viral load and label simultaneously ($\Sigma = 1$ for the processes $(S\rightarrow I_1)$, $(I_1\rightarrow I_2)$, $(I_2\rightarrow R)$, $(H_1\rightarrow H_2)$, $(H_2\rightarrow R)$) and $\Psi = 1$ if the label of the agents in compartment $X_t$ changes independently of the viral load (for birth and death processes and for the processes $(I_1\rightarrow H_1)$ and $(I_2\rightarrow H_2)$). $J_{X_t}$ is the new label of an agent performing a label switch independently of the viral load and previously labelled with $X_t$. Moreover, $V'_{X_t}$, $J^v_{X_t}$ are the new viral load and label of an agent with previous state $(X_t, V_t)$ in the case of a process in which the viral load and the label change simultaneously. Furthermore, we assume that $\Theta$ and $\Xi$ are two independent Bernoulli random variables describing whether a process happens ($\Theta = 1, \, \Xi = 1$) or not ($\Theta = 0, \, \Xi = 0$). We suppose that $P(\Theta = 1) = \lambda_{X_t} \Delta t$, being $\lambda_{X_t}$ the frequency of the microscopic process that rules the change of the microscopic variable $v$ for individuals labelled with $X_t$, while $P(\Xi = 1) = \lambda_{J_{X_t}, X_t}\Delta t$, where $\lambda_{J_{X_t}, X_t}$ is the frequency of the transition that causes the independent label switch from $X_t$ to $J_{X_t}$. In order for $P$ to be well defined, we must have that $\lambda_{X_t}\Delta t, \, \lambda_{J_{X_t}, X_t}\Delta t \le1$. The latter models the assumption according to which the larger the time interval, the higher the probability of having a label switch and/or a change of the viral load.

Independent label switch. Let us denote with $P(j\rightarrow i) = P(J_{X_t} = i|X_t = j)$ the conditional probability of switching from compartment $j$ to compartment $i$, with $i, \, j\in\mathcal{X}$, independently of a change of the microscopic state $v$. This probability concerns birth and death processes, and the isolation of infectious individuals, i.e. the label switches $(I_1\rightarrow H_1)$, $(I_2\rightarrow H_2)$. Specifically, we consider the following non–zero values for $P$:

● $P(B\rightarrow S) = {b}/{\rho_B(t)}\in[0, 1]$, where $b$ is a non–negative constant;

● $P(S \rightarrow D_\mu) = P(I_1\rightarrow D_\mu) = P(I_2\rightarrow D_\mu) = P(H_1\rightarrow D_\mu) = P(H_2\rightarrow D_\mu) = P(R\rightarrow D_\mu) = \mu\in[0, 1]$;

● $P(H_1 \rightarrow D_d) = P(H_2 \rightarrow D_d) = d\in[0, 1]$;

● $P(I_1 \rightarrow H_1) = P(I_2 \rightarrow H_2) = \alpha_H(v) \in [0, 1]$, where $\alpha_H(v)$ is an increasing function of the viral load $v$. It accounts for the fact that infectious people with a higher viral load are more likely to be identified. Indeed, performances of screening and diagnostic tests increase with the actual number of viral particles in the organism (see e.g. the interim guidance [32] on diagnostic testing for SARS–CoV–2). Moreover, for some infectious diseases a higher viral load is positively associated with a worse outcome and symptomatology (like for seasonal flu [21]).

We remark that, in principle, the probability of dying from the disease, $d$, may depend on the viral load $v$ like the probability of being isolated, $\alpha_H$. In the present work we assume that $d$ is constant because we are mainly interested in investigating the role of the viral load with respect to the isolation process.

The frequencies of the Markovian processes describing the switch between the different compartments may, in general, depend on both the departure and the arrival classes. It means that the process of switching from class $j$ to class $i$, that happens with probability $P(J_{X_t} = i|X_t = j)$, has frequency $\lambda_{J_{X_t}, X_t} = \lambda_{i, j}$. In particular, we consider:

● $\lambda_{S, B} = \lambda_b$, that is the frequency of new births or immigration;

● $\lambda_{D_\mu, j} = \lambda_\mu, \, j\in \mathcal{X}\setminus \{B, D_\mu, D_d\}$, that is the frequency of natural deaths;

● $\lambda_{D_d, H_1} = \lambda_{D_d, H_2} = \lambda_d$, that is the frequency of disease–induced deaths;

● $\lambda_{H_1, I_1}$, $\lambda_{H_2, I_2}$, that are the frequencies at which infectious individuals are isolated.

Simultaneous label switch. In our multi–agent system the first microscopic process causing simultaneously both a label switch and a progression of the viral load is the transition from susceptible ($S$) to infectious ($I_1$) state. This process has frequency $\lambda_{\beta}>0$. We express the corresponding transition probability as

that is the probability for an agent labelled with $X_t = S$ to change his/her label and zero viral load into $(I_1, v')$. Since this happens if a susceptible individual meets an infectious individual, we may regard $P(J_{X_t}^v = I_1, V'_{X_t} = v'|X_t = S)$ as a probability density distribution of the joint random variables $J_{X_t}^v$ and $V'_{X_t}$, given the probability $\rho_{I_1}(t)+\rho_{I_2}(t)$ of encountering an infectious individual, i.e.

that can be rewritten as

where $P(J^v_{X_t} = I_1|V'_{X_t} = v')$ is the probability density distribution of having an agent labelled with $I_1$ given that he/she has a viral load $v'$. In particular, $P(J^v_{X_t} = I_1|V'_{X_t} = v') = 1$ if $v'>0$, and $P(J^v_{X_t} = I_1|V'_{X_t} = v') = 0$ if $v' = 0$. $P(V'_{X_t} = v')$ is the probability density distribution of the random variable $V'_{X_t} = T_{\nu_\beta} v_0$ and it takes into account the microscopic rule describing the change of the state $v$ in terms of transition probabilities (see [23] for more details). It may be also expressed as $P(V'_{X_t} = v') = \nu_\beta P_S(v')$ where $P_S = \delta(v'-v_0)$.

Analogously, we may express the transition probabilities concerning the autonomous process and label switch as

where $P(J^v_{X_t} = i|V'_{X_t} = v')$ is the probability density distribution of having an agent labelled with $i$ given that he/she has a viral load $v'$, while $P_j(v\rightarrow v')$ is the transition probability describing the autonomous process of the viral load $v'$, given the previous viral load $v$, for agents labelled with $j$.

Remark 1. If $v'$ is such that $P(J^v_{X_t} = i|V'_{X_t} = v') = 1$ and $i = j$, then $P(J^v_{X_t} = i, V'_{X_t} = v'|X_t = j, V_t = v) = P_j(v\rightarrow v')$ is the transition probability that describes the change of the microscopic state $v'$ alone according to the rules (3)–(4) (see [23]).

In our case, the transitions to take into account are:

● $P(J_{X_t}^v = I_2, V'_{X_t} = v'|X_t = I_1, V_t = v) = P(J_{X_t}^v = I_2|V'_{X_t} = v')P_{I_1}(v\rightarrow v')$ and $P(J_{X_t}^v = H_2, V'_{X_t} = v'|X_t = H_1, V_t = v) = P(J_{X_t}^v = H_2|V'_{X_t} = v')P_{H_1}(v\rightarrow v')$, where $P(J_{X_t}^v = I_2|V'_{X_t} = v') = P(J_{X_t}^v = H_2|V'_{X_t} = v') = \eta(v')$. In principle, the probability $\eta(v')$ should increase by increasing the viral load $v'$, since individuals with a higher viral load are more likely to have reached the peak value. Here, for mathematical convenience, we approximate $\eta$ to the factor of viral load increase: $\eta = \nu_1\in[0, 1]$. Both $P_{I_1}(v\rightarrow v')$ and $P_{H_1}(v\rightarrow v')$ have average $v+\nu_1(1-v)$. In particular, we choose $P_{I_1}(v\rightarrow v') = P_{H_1}(v\rightarrow v') = \delta \big( v'- (v+\nu_1(1-v))\big)$.

● $P(J_{X_t}^v = R, V'_{X_t} = v'|X_t = I_2, V_t = v) = P(J_{X_t}^v = R|V'_{X_t} = v')P_{I_2}(v\rightarrow v')$ and $P(J_{X_t}^v = R, V'_{X_t} = v'|X_t = H_2, V_t = v) = P(J_{X_t}^v = R|V'_{X_t} = v')P_{H_2}(v\rightarrow v')$, where $P(J_{X_t}^v = R|V'_{X_t} = v') = \gamma(v')$ describes the probability for an infected individual of recovering. In principle, the probability $\gamma(v')$ should increase by decreasing the viral load $v'$, since individuals with lower viral load are more likely to have passed the infectious period. Similarly to what done for $\eta(v')$, we approximate this probability to the factor of viral load decay: $\gamma = \nu_2\in[0, 1]$. $P_{I_2}(v\rightarrow v')$ and $P_{H_2}(v\rightarrow v')$ have average $v-\nu_2 v$. In particular, we choose $P_{I_2}(v\rightarrow v') = P_{H_2}(v\rightarrow v') = \delta\big(v'-v(1-\nu_2)\big)$.

● $P(J_{X_t}^v = I_1, V'_{X_t} = v'|X_t = I_1, V_t = v) = P_{I_1}(v\rightarrow v')$, $P(J_{X_t}^v = H_1, V'_{X_t} = v'|X_t = H_1, V_t = v) = P_{H_1}(v\rightarrow v')$, while $P(J_{X_t}^v = I_2, V'_{X_t} = v'|X_t = I_2, V_t = v) = P_{I_2}(v\rightarrow v')$, $P(J_{X_t}^v = H_2, V'_{X_t} = v'|X_t = H_2, V_t = v) = P_{H_2}(v\rightarrow v')$.

● $P(J_{X_t}^v = R, V'_{X_t} = v'|X_t = R, V_t = v) = P_{R}(v\rightarrow v') = \delta\big(v'-v(1-\nu_2)\big)$.

The frequency of these transitions is the frequency of the corresponding microscopic process, i.e. $\lambda_{X_t} = \lambda_{\gamma}$ for $X_t \in \lbrace I_1, I_2, H_1, H_2, R \rbrace$.

3.   Aggregate description: from kinetic to macroscopic equations

The kinetic equations describing the evolution of $f_i(t, v), \, i \in \mathcal{X}$, can be derived in the same way as in [24]. Namely, the system of the weak equations for the $f_i$'s is the following:

$i\in\mathcal{X}$, where $\varphi:[0, 1]\rightarrow \mathbb{R}$ is a test function. In (6), the second and third lines account for the Markovian processes describing the label switches that happen, respectively, independently of the evolution of the viral load, and simultaneously with the evolution of the viral load. The frequency $\lambda_{i}$ of the Markovian process due to an interaction with a background corresponds to the frequency of changing the microscopic state $v$ and it is, as previously stated, $\lambda_i = \lambda_\gamma$, $\forall i \in \lbrace I_1, I_2, H_1, H_2, R \rbrace$.

From (6), we derive the kinetic equations describing the evolution of the distribution functions $f_i$'s, $i\in\mathcal{X}$. For $i\in\mathcal{X}\setminus\{B, D_d, D_\mu\}$, namely the classes of living individuals, we get:

● susceptible individuals ($i = S$)

● infectious individuals with increasing viral load ($i = I_1$)

● infectious individuals with decreasing viral load ($i = I_2$)

● isolated individuals with increasing viral load ($i = H_1$)

● isolated individuals with decreasing viral load ($i = H_2$)

● recovered individuals ($i = R$)

Equations (7)–(12) have to hold for every $\varphi: [0, 1]\to\mathbb{R}$.

In order to obtain the equations for the macroscopic densities and viral load momentum of each compartment, we set $\varphi(v) = v^n$ in (7)–(12), with $n = 0, 1$, respectively. Since setting $\varphi(v) = v^n$ in the evolution equations (7)–(12) leads to the appearance of the $(n+1)$th moment of $f_i$, namely $\int_0^1f_i(t, v)v^{n+1} dv$, we need to find a closure. Specifically, for each compartment we consider a monokinetic closure in the form

i.e. we assume that all the agents of the same compartment at a given time $t$ have the same viral load.

Remark 2. We consider that a monokinetic closure is appropriate in this framework, as the microscopic rules describing the evolution of the viral load (3)–(4) are deterministic, i.e. there is no diffusion related to stochastic fluctuations.

As already observed, if $\rho_i(t) = 0$, then the mean viral load is not well defined. Notwithstanding, since $f_i$ is defined as a Dirac delta, we have that if $\varphi(v)$ is a test function, then

Then, we consider test functions such that

When $\varphi(v) = 1$ or $\varphi(v) = v$, namely the test functions allowing to recover the densities and momentum, respectively, condition (14) is satisfied. Moreover, we deal with terms in the form

with $\psi = \lambda_\gamma\eta$, $\psi = \lambda_\gamma\gamma$ and $\psi = \lambda_{H_j, I_j}\alpha_H$, $j = 1, 2$. Since we assumed that the probabilities $\eta, \, \gamma$ are constant, then the integral (15) with $\psi = \lambda_\gamma\eta$ or $\psi = \lambda_\gamma\gamma$ is well defined – i.e. it is not divided by a vanishing density – for both test functions $\varphi(v) = 1$ and $\varphi(v) = v$. As far as the isolation terms are concerned, i.e. $\psi(v) = \lambda_{H_j, I_j}(t)\alpha_H(v)$, $j = 1, 2$, we have that, applying the monokinetic closure, the integral (15) reads

Hence, we have to choose the isolation frequency and probability function in such a way that the latter quantity is well defined.

The macroscopic model is given by the following system of non–linear ordinary differential equations:

For convenience of notation, in (16) we have denoted with the upper dot the time derivative and omitted the explicit dependence on time of the state variables.

To model (16) we associate the following generic initial conditions

Remark 3. We observe that, by assuming $\gamma(v) = \nu_2$ and $\eta(v) = \nu_1$ (see Section 2.3), we are not keeping into account the dependence of the transitions $(I_2\rightarrow R)$, $(H_2\rightarrow R)$, $(I_1\rightarrow I_2)$, $(H_1\rightarrow H_2)$ on the viral load value. However, with these choices, the recovery rate of the infected individuals in classes $I_2$ and $H_2$ is $\lambda_{\gamma}\nu_2$, that is the decay rate of the viral load. Analogously, the rate of transition from $I_1$ [resp. $H_1$] to $I_2$ [resp. $H_2$] is the increase rate of the viral load, $\lambda_{\gamma}\nu_1$.

As far as the products $\lambda_{H_j, I_j}(t)\alpha_H(v)$, $j = 1, 2$, are concerned, let us note that, if both the frequencies $\lambda_{H_j, I_j}$, $j = 1, 2$, and the probability $\alpha_H(v)$ are assumed to be constant, then from (16) we retrieve a classical SIR–like model with constant isolation rate. The qualitative analysis of the ensuing model can be easily obtained and is here omitted.

We focus instead on the impact of viral load–sensitivity of tests and frequency of testing activities on the epidemic dynamics and consider the case that:

● the probability for infectious individuals to be isolated, $\alpha_H(v)$, linearly increases with their viral load: $\alpha_H(v) = \alpha v$, where $\alpha\in[0, 1]$ is a constant;

● the frequencies $\lambda_{H_j, I_j}(t)$ are linearly dependent on the densities of infectious individuals: $\lambda_{H_j, I_j}(t) = \lambda_\alpha\rho_{I_j}(t)$, $j = 1, 2$, where $\lambda_\alpha$ is a non–negative constant. Namely, we assume that the efforts made by public health authorities in screening and diagnostic activities increase with the increase of infectious presence in the community. Indeed, when infectious individuals are few, search activities could be highly expensive and little effective.

With these choices, in system (16), the isolation terms become

Equilibria and stability properties of model (16)–(18) will be investigated in the following section.

4.   Qualitative analysis

Since in model (16)–(18) the differential equations for $\rho_S,$ $\rho_{I_1},$ $\rho_{I_2},$ $n_{I_1},$ $n_{I_2}$ do not depend on $\rho_{H_1},$ $\rho_{H_2},$ $\rho_{R},$ $n_{H_1},$ $n_{H_2},$ $n_{R}$, it is not restrictive to limit our analysis to system

It is straightforward to verify that the region

with initial conditions in (17) is positively invariant for model (19), namely any solution of (19) starting in $\mathcal{D}$ remains in $\mathcal{D}$ for all $t\geq 0$.

In the following, we search for model equilibria and derive suitable thresholds ruling their local or global stability.

4.1.   Disease–free equilibrium and its stability

The model (19) has a unique disease–free equilibrium (DFE), given by

It is obtained by setting the r.h.s. of equations (19) to zero and considering the case $\rho_{I_1} = \rho_{I_2} = 0$.

Proposition 1. The DFE of system (19) is locally asymptotically stable (LAS) if $\mathcal{R}_0<1$, where

Otherwise, if $\mathcal{R}_0>1$, then it is unstable.

Proof. The Jacobian matrix $J$ of system (19) evaluated at the DFE reads

with

One can immediately get the eigenvalues $l_1 = -\lambda_\mu \mu<0$, $l_2 = J_4<0$, $l_3 = J_{5}<0$, while the other two are determined by the submatrix

From the sign of the entries of $\bar J$, it follows that det$(\bar J)\geq0$ implies tr$(\bar J)<0$. Hence, if det$(\bar J)>0$ or, equivalently, if $\mathcal{R}_0<1$, with $\mathcal{R}_0$ given in (21), then the DFE is LAS. Otherwise, if $\mathcal{R}_0>1$, then it is unstable.

The threshold quantity $\mathcal{R}_0$ is the so–called basic reproduction number for model (19), 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. One can easily verify that the same quantity can be obtained as the spectral radius of the so–called next generation matrix [29].

As far as the global stability of the DFE, we prove the following theorem.

Theorem 4.1. The DFE of system (19) is globally asymptotically stable if $\mathcal{R}_0< 1$.

Proof. Consider the following function

It is easily seen that the $\mathcal{L}$ is non–negative in $\mathcal{D}$ (see (20)) and also $\mathcal{L} = 0$ if and only if $\rho_{I_1} = \rho_{I_2} = 0$. The time derivative of $\mathcal{L}$ along the solutions of system (19) in $\mathcal{D}$ reads

It follows that $\dot{\mathcal{L}}\leq0$ for $\mathcal{R}_0 < 1$ with $\dot{\mathcal{L}} = 0$ only if $\rho_{I_1} = \rho_{I_2} = 0$. Hence, ${\mathcal{L}}$ is a Lyapunov function on $\mathcal{D}$ and the largest compact invariant set in $\{\left(\rho_S, \rho_{I_1}, \rho_{I_2}, n_{I_1}, n_{I_2}\right)\in \mathcal{D}: \dot{\mathcal{L}} = 0 \}$ is the singleton {$DFE$}. Therefore, from the La Salle's invariance principle [19], every solution to system (19) with initial conditions in (17) approaches the DFE, as $t\rightarrow +\infty$.

As an alternative proof, one may adopt the approach developed by Castillo–Chavez et al. in [3].

4.2.   Endemic equilibria

Let us denote with

the generic endemic equilibrium of model (19), obtained by setting the r.h.s. of equations (19) to zero and considering the case $\rho_{I_1}+\rho_{I_2}>0$. Note that, if it were $\rho_{I_1}^E = 0$ [resp. $\rho_{I_2}^E = 0$], then from (19c) it would follow that $\rho_{I_2}^E = 0$ [resp. $\rho_{I_1}^E = 0$]. Hence, it must be $\rho_{I_1}^E, \, \rho_{I_2}^E>0$.

More precisely, by rearranging equations (19a)–(19b)–(19c)–(19d), one obtains

Substituting the expressions (22) into (19e), one gets $n_{I_1}^E$ as a positive root of the equation

Due to the complexity of equation (23), we renounce to get an explicit expression for $n_1^E$ and, hence, to derive the existence conditions and number of endemic equilibria. However, we will make use of bifurcation analysis to show that a unique branch corresponding to a unique endemic equilibrium emerges from the criticality, namely at DFE and $\mathcal{R}_0 = 1$.

4.3.   Central manifold analysis

To derive a sufficient condition for the occurrence of a transcritical bifurcation at $\mathcal{R}_0 = 1$, we can use a bifurcation theory approach. We adopt the approach developed in [10,29], which is based on the general center manifold theory [16]. 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:

where

and

Note that in (24) and (25) the product $\lambda_\beta\nu_\beta$ has been chosen as bifurcation parameter, $\overline{\lambda_\beta\nu_\beta}$ is the critical value of $\lambda_\beta\nu_\beta$, $\mathbf{x} = \left(\rho_S, \rho_{I_1}, \rho_{I_2}, n_{I_1}, n_{I_2}\right)$ is the state variables vector, $\mathbf{F}$ is the right–hand side of system (19), and $\mathbf{z}$ and $\mathbf{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 $\lambda_\beta\nu_\beta = \overline{\lambda_\beta\nu_\beta}$).

Observe that $\mathcal{R}_0 = 1$ is equivalent to:

so that the disease–free equilibrium is stable if $\lambda_\beta\nu_\beta<\overline{\lambda_\beta\nu_\beta}$, and it is unstable when $\lambda_\beta\nu_\beta>\overline{\lambda_\beta\nu_\beta}$.

The direction of the bifurcation occurring at $\lambda_\beta\nu_\beta = \overline{\lambda_\beta\nu_\beta}$ can be derived from the sign of coefficients (24) and (25). More precisely, if $A>0$ [resp. $A<0$] and $B>0$, then at $\lambda_\beta\nu_\beta = \overline{\lambda_\beta\nu_\beta}$ there is a backward [resp. forward] bifurcation.

For our model, we prove the following theorem.

Theorem 4.2. System (19) exhibits a forward bifurcation at DFE and $\mathcal{R}_0 = 1$.

Proof. From the proof of Proposition 1, one can verify that, when $\lambda_\beta\nu_\beta = \overline{\lambda_\beta\nu_\beta}$ (or, equivalently, when $\mathcal{R}_0 = 1$), the Jacobian matrix $J(DFE)$ 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 $\mathbf{z\cdot}\mathbf{w} = 1$ are:

with

and

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 $\mathbf{F}$, it follows that:

and

where $z_2, \, z_3, \, w_2, \, w_3, \, w_4, \, w_5>0$ and $w_1<0$. Then, $A<0<B$. Namely, when $\lambda_\beta\nu_\beta-\overline{\lambda_\beta\nu_\beta}$ changes from negative to positive, DFE changes its stability from stable to unstable; correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable. This completes the proof.

5.   Numerical simulations

In this section, we present and compare some numerical solutions of both the stochastic particle model (5) and the macroscopic model (16).

Our aim is to qualitatively assess the interplay between the evolution of individuals' viral load and the disease spread and isolation control. Hence, demographic and epidemiological parameters values do not address a specific infectious disease and/or spatial area. They refer to a generic epidemic outbreak where control strategies rely on isolation of infectious individuals, as typically happens for new emerging infectious diseases (e.g., 2003–2004 SARS outbreak [31], 2014–2016 Western African Ebola virus epidemic [4], the first phase of the ongoing COVID–19 pandemic [33]).

Numerical simulations are performed in Matlab^Ⓡ [26]. We implement a Monte Carlo algorithm to simulate the stochastic particle model (5) and the 4th order Runge–Kutta method with constant step size for integrating the system (16). Platform–integrated functions are used for getting the plots.

5.1.   Parametrization

The time span of our numerical simulations is set to $t_f = 1$ year. We are considering an SIR–like model with demography and constant net inflow of susceptibles $\lambda_b b$. Since travel restrictions are usually implemented during epidemic outbreaks, we assume that $\lambda_b 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, $\bar N$ denotes the total resident population at the beginning of the epidemic, and $N_{tot}$ is the total system size. Note that $N_{tot}$ accounts for agents belonging to all model compartments $\mathcal{X}$ (including $B$, $D_\mu$, $D_d$), whereas $\bar N$ refers only to living individuals.

We assume a population of $\bar N = 10^6$ individuals, representing, for example, the inhabitants of a European metropolis. Fluctuations in a time window of just over a year are considered negligible. The most recent data by European Statistics refer to 2019 and provide an average crude birth rate $b_r = 9.5/1,000$ years$^{-1}$ [13] and an average crude death rate $\lambda_\mu\mu = 10.2/1,000$ years$^{-1}$ [12]. The total (constant) system size $N_{tot}$ is set to $N_{tot} = \bar N/(1-\lambda_bbt_f)$, in such a way $N_{tot} = \bar N +\lambda_bbt_f N_{tot}$ is given by the sum of the initial population, $\bar N$, and the total inflow of individuals during the time span considered, $\lambda_bbt_f N_{tot}$.

For the epidemiological parameters we take the following baseline values:

In particular, the product $\lambda_\gamma\nu_1$ can be interpreted as the inverse of the average time from exposure to viral load peak, whilst $\lambda_\gamma\nu_2$ as the inverse of the average time from viral load peak to recovery. The disease–induced death rate $\lambda_d d$ is estimated through the formula given by Day [7]:

where $C_F$ is the fatality rate and $T$ is the expected time from isolation until death. We assume $C_F = 1$% and $T = 1/(\lambda_\gamma\nu_2) = 10$ days. As far as the initial viral load of infected individuals, $v_0$, is concerned, we assume that it is 1% of the maximum reachable value ($v = 1$), namely $v_0 = 0.01$. Finally, for the Monte Carlo simulation of the particle model (5), we further assume $\lambda_b = \lambda_{\beta} = 1$ days$^{-1}$ and $\lambda_\mu = \lambda_d = 0.01$ days$^{-1}$.

Initial data are set to the beginning of the epidemic, namely we consider a single infectious individual in a totally susceptible population:

All the parameters of the model as well as their baseline values are reported in Table 1.

5.2.   The uncontrolled epidemic outbreak

First, we numerically investigate the impact of the epidemiological parameters on the basic reproduction number $\mathcal{R}_0$, see (21). By considering the baseline parameters values, we obtain that the ratio between $\mathcal{R}_0$ and the transmission rate $\lambda_\beta\nu_\beta$ is about 13.83. Fig. 1 displays the contour plot of $\mathcal{R}_0$ versus $\lambda_\gamma\nu_1$ (x–axis values) and $\lambda_\gamma\nu_2$ (y–axis values). We vary the average period of viral load increase, $1/(\lambda_\gamma\nu_1)$, in the range $[2, 14]$ days and the average period of viral load decay, $1/(\lambda_\gamma\nu_2)$, in the range $[5, 25]$ days. We obtain that $\mathcal{R}_0$ decreases with both $\lambda_\gamma\nu_1$ and $\lambda_\gamma\nu_2$, from a maximum of $\mathcal{R}_0 = 10.39$ for $\lambda_\gamma\nu_1 = 1/14$ days$^{-1}$ and $\lambda_\gamma\nu_1 = 1/25$ days$^{-1}$ to a minimum of $\mathcal{R}_0 = 1.87$ for $\lambda_\gamma\nu_1 = 1/2$ days$^{-1}$ and $\lambda_\gamma\nu_2 = 1/5$ days$^{-1}$.

Let us now set the basic reproduction number to the baseline value $\mathcal{R}_0 = 4$ and investigate the epidemic dynamics in absence of isolation control ($\alpha_H\equiv 0$). Numerical simulations are displayed in Fig. 2. We compare densities and viral load means. Specifically, solid lines refer to the solutions of the macroscopic model (16) and markers to those of the stochastic particle model (5). We note a good match between the two approaches in predicting the dynamics of compartment sizes $\rho_iN_{tot}$, $i\in\mathcal{X}$ (grey scale colour): an epidemic outbreak invades the population, by reaching a prevalence peak of approximately $(\rho_{I_1}+\rho_{I_2})N_{tot} = 531,000$ in 61 days; after 1 year the prevalence is almost zero and susceptible individuals are just about $21,700$. On the contrary, the dynamics of compartment mean viral loads $n_i/\rho_i$, $i\in\mathcal{X}$ (blue scale colour), may be different in the particle w.r.t. the macroscopic model: the match is good as long as the corresponding compartment size is not so small to make the effect of stochasticity relevant. This discrepancy is evident at the beginning [resp. at the end] of the time horizon, namely when the density of infectious and recovered [resp. of infectious] individuals is almost zero, see Fig. 2(b)–(c)–(d) [resp. Fig. 2(c)–(d))]. In particular, from Fig. 2(c)–(d), we see that, at the end of the epidemic wave, according to the particle model (blue markers) the mean viral loads of infectious individuals fluctuate until approaching zero when the corresponding compartment becomes empty. Instead, the macroscopic model predicts that the same means remain approximately constant at a positive value (blue solid lines), suggesting that the first moment $n_{I_1}$ [resp. $n_{I_2}$] and the density $\rho_{I_1}$ [resp. $\rho_{I_2}$] go to zero with the same speed. This is due to the inconsistency of average quantities, like the mean viral loads, when the number of particles is very small. In that case, the deterministic macroscopic model cannot be justified by means of the law of large numbers and statistical fluctuations must be taken into account.

5.3.   Viral load–dependent vs. constant isolation control

Here, we introduce the isolation control in the epidemic model and assess how the frequency of testing and the viral load–sensitivity of tests can affect epidemic dynamics. To this aim, we compare the following simulation scenarios:

S1 viral load–dependent isolation, as studied here: $\lambda_{H_j, I_j}(t) = \lambda_\alpha\rho_{I_j}(t)$, $j = 1, 2$, and $\alpha_H(v) = \alpha v$;

S2 constant isolation, as in classical epidemic models: $\lambda_{H_j, I_j}(t) = \lambda_\alpha$, $j = 1, 2$, and $\alpha_H(v) = \alpha$.

In order to make the two scenarios properly comparable, we make the following considerations. In the case S2, the product $\lambda_\alpha\alpha$ represents the rate at which infectious individuals are isolated in the unit of time. In the case S1, in the microscopic model, the same rate is given by $\lambda_\alpha\alpha$ multiplied by the individual microscopic viral load $v$ and by the densities $\rho_{I_j}$, $j = 1, 2$; whereas, in the macroscopic model (16), due to the monokinetic closure, this rate is given by $\lambda_\alpha\alpha$ multiplied by the momentum $n_{I_j}$, $j = 1, 2$, related to the corresponding compartment (see (18)). Thus, we assume that the value of $\lambda_\alpha\alpha$ in scenario S1 (say, $\left.\lambda_\alpha\alpha\right.|_{S1}$) is given by the value adopted in scenario S2 ($\left.\lambda_\alpha\alpha\right.|_{S2}$) rescaled by a normalization factor $M$:

where $M$ represents an average quantity for the $n_{I_j}$'s, $j = 1, 2$. In order to estimate $M$, we consider the model (16) in absence of isolation control ($\alpha_H\equiv 0$) and denote by $n_{I_1}^{unc}(t)$ and $n_{I_2}^{unc}(t)$ the corresponding solutions for $n_{I_1}(t)$ and $n_{I_2}(t)$, respectively. Then, $M$ is set to

where $\bar{n}_{I_j}$ are the average values of $n_{I_j}^{unc}(t)$ over $[0, t_f]$, namely

The numerical value for $\left.\lambda_\alpha\alpha\right.|_{S2}$ is set to $\left.\lambda_\alpha\alpha\right.|_{S2} = 0.1$ days$^{-1}$. For the Monte Carlo simulation we set $\lambda_\alpha = 15$ days$^{-1}$ in scenario S1 and $\lambda_\alpha = 1$ days$^{-1}$ in scenario S2.

Figs. 3 and 4 display the numerical simulations in scenarios S1 (grey scale colour) and S2 (blue scale colour). Specifically, solid lines refer to the solutions of the macroscopic model (16) and markers to those of the stochastic particle model (5). As far as the match between the two approaches is concerned, we note that in the case of constant control S2 considerations similar to those made in Section 5.2 apply. Instead, in the case of viral load–dependent control S1, solutions by particle and macroscopic models are qualitatively similar but quantitatively different. This is an expected result because the derivation of the macroscopic model relies on an approximation through the monokinetic closure (13), which acts by levelling the viral loads of all agents belonging to a given class to their average value. However, notwithstanding the postulated monokinetic closure, the matching is quite good, as the peak given by the macroscopic model is only mildly underestimated. Again, we remark that the macroscopic model may not well reproduce the compartment mean viral loads as predicted by the stochastic model (Fig. 4). This happens at the beginning and at the end of the time horizon when some compartments are almost empty (e.g., those of infectious and isolated individuals). In such cases, the law of large numbers does not apply and the concept of theoretical mean departs consistently from that of empirical mean.

As far as the comparison between scenarios S1 and S2 is concerned, from Fig. 3 we note that in the first case (grey scale colour) the epidemic outbreak occurs earlier and with a lower peak w.r.t. the second case (blue scale colour), but the tails of the infected curves are longer. In order to investigate these differences more deeply, we consider the solutions by the macroscopic model (16) and report in Table 2 some relevant epidemiological quantities, including the value of infectious prevalence peak and the time it occurs, and the endemic value of infectious prevalence, $(\rho_{I_1}^E+\rho_{I_2}^E)N_{tot}$. In scenario S1, the endemic components $\left.\rho_{I_1}^E\right|_{S1}$, $\left.\rho_{I_2}^E\right|_{S1}$ are computed through the expressions in (22). In particular, in our numerical set, equation (23) admits three positive roots, but just one of them makes all the other endemic components (22) positive, hence a unique endemic equilibrium exists. In scenario S2, one can easily verify that the unique endemic equilibrium has

We also compute the value at the final time $t_f = 1$ year of three cumulative quantities: the cumulative incidence CI$(t)$, i.e. the total number of new cases in $[0, t]$; the cumulative isolated individuals, i.e. the total number of infectious individuals that tested positive in $[0, t]$, and the cumulative deaths CD$(t)$, i.e. the disease–induced deaths in $[0, t]$. In our model we have, respectively:

From Table 2, we note that the epidemic peak in scenario S1 is almost halved compared to the scenario S2 and occurs 36 days before. By contrast, the endemic infectious prevalence is much greater in scenario S1 w.r.t. S2: 289 vs. 76. Interestingly, the differences in the cumulative quantities CI$(t_f)$, CH$(t_f)$, CD$(t_f)$ are, instead, minimal: in case of viral load–dependent isolation the cumulative incidence at 1 year is approximately 2% greater than the corresponding quantity in the case of constant isolation, while the cumulative isolated individuals [resp. deaths] are about 0.2% greater [resp. smaller].

The viral load–dependent isolation function reflects the assumption that an infectious individual with high viral load is more likely to be identified: it may represent the efficiency of the test that, according to its sensitivity, is capable of detecting different concentrations of virus particles per $ml$ [20]. Assuming a constant isolation function means, instead, that all infectious individuals have the same probability of being detected and diagnosed. In other words, infectious individuals with sufficiently low [resp. high] viral load have a probability of being diagnosed higher [resp. lower] in scenario S2 with respect to S1.

5.4.   Tracking individuals' viral load

One of the advantages of a particle model is the possibility to track the trends of all the agents of the system. Here, we are interested in tracking the evolution of individuals' viral load during the simulation time span. To this aim, we consider the particle model (5) with viral load–dependent isolation control (scenario S1) and retrieve the viral load evolution of every single agent. In Fig. 5, we report the temporal dynamics of $v$ for five selected agents, who show different courses of the disease. Different line markers and/or colours refer to the different epidemiological compartments the agents pass through; the meaning is specified in the figure legend. Note that two agents die after having acquired the infection: one of natural causes (first curve from the left), the other one from the disease (second curve from the left). The other three agents survive and finally recover from the infection: two of them are identified and isolated during the infectious period (third and fourth curve from the left), while the last one remains free to move (fifth curve from the left). We also remark that individuals may recover before their viral load becomes null and that it may take a long time after recovering for $v$ to completely vanish. Such a trend is linked to the choice of a constant probability of recovery, $\gamma = \nu_2$. This is in agreement with experimental observations of viral load curves, that show that individuals are no longer infectious before the complete disappearance of the virus [5,6,15,17,21]. Nonetheless, the mathematical description could be refined by setting a $v-$dependent and decreasing probability of recovering, $\gamma(v)$.

6.   Discussion and conclusion

In this work, we have proposed a microscopic stochastic model allowing one to describe the spread of an infectious disease through social contacts. Each individual is identified by the epidemiological compartment to which he/she belongs and by his/her viral load. Binary interactions between susceptible and infectious individuals may cause the susceptible to acquire a positive viral load $v$ and, as a consequence, to get infected. The viral load progression due to physiological processes is different from person to person, it determines the health status of the individuals and, therefore, the epidemiological compartment to which they belong. In particular, we have here considered the case that the viral load influences explicitly the isolation mechanism, i.e. the switch from $I_i$ to $H_i$, $i = 1, 2$. In this sense, the present work manages to deal with the heterogeneity of the individuals' viral load, that is explicitly encoded in the individual viral load progression and in the probability of being diagnosed.

We have derived from the stochastic particle model a kinetic description by means of evolution equations of the distribution of the viral load in each compartment. Finally, by making use of a monokinetic closure, we have obtained a macroscopic description. The ensuing macroscopic model is a system of non–linear ordinary differential equations for the macroscopic densities and viral load momentum of the compartments. We have performed a qualitative analysis allowing to state that our system has a unique disease–free equilibrium (DFE) that is globally asymptotically stable if $\mathcal{R}_0<1$ and that the system (16) exhibits a transcritical forward bifurcation at DFE and $\mathcal{R}_0 = 1$.

Our numerical tests have allowed us to compare the predictions yielded by the particle and the macroscopic models. In particular, concerning the densities of the compartments, the numerical tests have confirmed the matching between the two approaches in the case that the probability of being isolated, $\alpha_H$, is constant, thereby validating the macroscopic model as a reliable approximation of the particle model more amenable to analytical investigations and quick numerical solutions. In the case of a viral load–dependent isolation, i.e. $\alpha_H(v) = \alpha v$, along with a density–dependent frequency of testing, we have seen that the qualitative trends of the compartment densities predicted by the particle (5) and the macroscopic (16) models are close but do not coincide exactly: this is a consequence of the approximation made through the monokinetic closure. Concerning instead the mean viral loads of the compartments, they are well reproduced when the density of the compartment is high enough: this is a consequence of the law of large numbers. Instead, when the compartment is almost empty and there are few incoming and outgoing individuals, the statistical average described by the macroscopic model is no longer reliable. Therefore, while the macroscopic model permits quicker numerical solutions, the particle model allows one to compute more accurately the viral load of the individuals. Moreover, the particle model allows one to track the viral load of every single agent and to investigate different possible individual evolutions of the disease.

Deliberately, we have not tried to match real scenarios by calibrating or comparing the results of our models with empirical data. In fact, our aim was first to propose a simple compartmental model including the viral load as microscopic variable. As a consequence, we wanted to explore prototypical scenarios and to compare them to those predicted by classical epidemic models, by focusing on the impact of having a viral load–dependent isolation in place of a constant isolation rate. We have seen that in the case of a viral load–dependent isolation the epidemic outbreak occurs earlier and with a lower peak (almost halved) w.r.t. to the constant isolation case. However, the cumulative disease–related quantities one year after the onset of the epidemic are comparable, while the endemic infectious prevalence is much greater in the viral load–dependent isolation scenario. This may be explained in terms of the viral load–sensitivity and frequency of the testing activities that are embodied in the choice of the functions $\alpha_H(v)$ and $\lambda_{H_1, I_1}(t)$, $\lambda_{H_2, I_2}(t)$, respectively.

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 (more amenable, of course, to analytical and numerical investigations), that provides for a richer description of the disease spreading in the host population. Here we only considered the explicit influence of the viral load on the isolation mechanism, but, in principle, other switches of individuals between compartments may depend on the viral load at the microscopic level, and on the viral load mean at the macroscopic level. Therefore, more complex situations, such as super–spreading events, that have been proved to be of the utmost importance for example during the COVID–19 pandemic [15], could be addressed. The heterogeneity of transmission could be included by making the disease transmission rate from infectious to susceptible individuals dependent on the viral load. Also, in such a way, different initial viral loads of the infectious individual first introduced in the community may give rise to different epidemic scenarios.

Acknowledgments

This work was supported by GNFM (Gruppo Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di Alta Matematica) through a 'Progetto Giovani 2020' grant. This work was also partially supported by the Italian Ministry of University and Research (MUR) through the "Dipartimenti di Eccellenza" Programme (2018-2022), Department of Mathematical Sciences "G. L. Lagrange", Politecnico di Torino (CUP: E11G18000350001). This work is also part of the activities of the PRIN 2017 project (No. 2017YBKNCE) "Multiscale phenomena in Continuum Mechanics: singular limits, off–equilibrium and transitions", and of the PRIN 2020 project (No. 2020JLWP23) "Integrated Mathematical Approaches to Socio–Epidemiological Dynamics".