A Discrete Model for the Evolution of Infection Prior to Symptom Onset

Abstract

We consider a between-host model for a single epidemic outbreak of an infectious disease. According to the progression of the disease, hosts are classified in regard to the pathogen load. Specifically, we are assuming four phases: non-infectious asymptomatic phase, infectious asymptomatic phase (key-feature of the model where individuals show up mild or no symptoms), infectious symptomatic phase and finally an immune phase. The system takes the form of a non-linear Markov chain in discrete time where linear transitions are based on geometric (main model) or negative-binomial (enhanced model) probability distributions. The whole system is reduced to a single non-linear renewal equation. Moreover, after linearization, at least two meaningful definitions of the basic reproduction number arise: firstly as the expected secondary asymptomatic cases produced by an asymptomatic primary case, and secondly as the expected number of symptomatic individuals that a symptomatic individual will produce. We study the evolution of infection transmission before and after symptom onset. Provided that individuals can develop symptoms and die from the disease, we take disease-induced mortality as a measure of virulence and it is assumed to be positively correlated with a weighted average transmission rate. According to our findings, transmission rate of the infection is always higher in the symptomatic phase yet under a suitable condition, most of the infections take place prior to symptom onset.

1. Introduction

Epidemiological models are a powerful mathematical tool to study the spread of infectious diseases, [1,2,3,4,5,6]. They can be used to predict the fraction of hosts in a population who eventually will become infected during each epidemic outbreak of a transmissible disease. Moreover, single-outbreak models can be concatenated to study a sequence of waves for a disease like COVID-19. The epidemic model introduced is intended for human or animal populations and for diseases with some sort of infection prior to symptom onset. Examples of diseases with asymptomatic carriers are typhoid, HIV, C. difficile, influenza, cholera, tuberculosis and COVID-19. Some people can be infected with flu viruses, for instance, and have no symptoms but may still be able to spread the virus to their close contacts. See also [7] for a paper on the emergence and spread of SARS-CoV-2 variants of concern.

We choose discrete time for the ease of the implementation of the models, [8,9]. It could be said that continuous-time and discrete-time models are the two faces of the same coin, the former being more popular in general. However, the leverage of the latter is that you do not need to discretize when doing numerical simulations and the probabilistic background of the models are easier to uncover. Discrete-time models are typically finite dimensional, as the one considered in here since the we assumed six different disease stages including two degrees of symptoms. However, discrete-time models can be also infinite dimensional when dealing with continuously structuring variables. See, e.g., [10] for the theory and implementation of Integral Projection Models, which are useful for studying ecological problems.

On the other hand, epidemiological models are the base to undertake evolutionary studies of a disease or a pathogen, [11,12,13,14]. From the epidemiological point of view, one could think that the appearance of symptoms of a disease enhances or depresses transmission success, [14,15,16,17]. On the one hand, symptoms could enhance transmission, for example if coughing and sneezing increases spread beyond what would be predicted for a given pathogen load. On the other hand, symptoms could alter host behavior to reduce transmission, e.g., via isolation or avoidance of individuals showing obvious signs of infection. Interestingly, both effects can happen together and we wonder which one is most important, that is, if one of them dominates the other. To answer this question we follow a simple approach consisting in the maximization of the fitness of the pathogen (e.g., basic reproduction number of the disease, [18]) under a tradeoff between disease-induced mortality (virulence, [19]) and transmission of the infection. More precisely, by evolution in the present paper we mean the competition of pathogens with different virulence levels and accordingly different infection transmission, as virulence is assumed correlated with the average transmission rate along infection phases, and all other traits and parameters of the model are fixed. Moreover, it is assumed that evolution leads to the maximization of the reproduction number of infectives rather than other epidemiological indicators like the initial growth rate of the infection. At evolutionary equilibrium, one would expect an intermediate modest level of virulence.

See for instance [12,13,19,20,21] for other evolutionary studies using tradeoffs. For evolutionary epidemiology using the adaptive dynamics approach see, e.g., [11,17]. It is worth to mention that to check the robustness of the outcomes, we have used three different definitions of the basic reproduction number, giving the same optimal virulence value, [19]. Those definitions depend on what is understood as infection event.

For a systematic approach to the analytical/numerical computation of the basic reproduction number in general structured population models, see [22,23] and the references therein.

2. The Model: Markov Chain

With the aim of studying the evolution of infectious processes with transmission in two different phases: a first phase with mild or no symptoms and a second phase where hosts are showing up symptoms, let us consider a discrete-time model for each epidemic outbreak. According to the viral load, we deal with hosts classified into the following phases: non-infectious asymptomatic phase, infectious asymptomatic phase (mild or no symptoms), infectious symptomatic phase and immune phase (natural immunity). For simplicity, no demographic turnover is considered and neither non-pharmaceutical measures nor vaccination are explicitly taken into account. Moreover, reinfections, that is, the loss of immunity during the epidemic outbreak, will be considered at the end of the paper.

The system takes the form of a non-linear model in discrete-time $t=0,1,2,\cdots$ days, as a Markov chain for the fraction of Susceptible $S_t$, Exposed $E_t$ (latent who are not infectious yet), infectious Asymptomatic $A_t$, Infectious symptomatic $I_t$, Removed $R_t$ (alive and immune) and (disease-related) Deceased $D_t$ hosts. See the flow diagram in Figure 1 for an explanation of the transitions considered in the present model. In words: susceptible hosts, $S_t$, enter into the exposed class, $E_t$, upon exposition to the pathogen by a contact with infectious hosts, $A_t$ or $I_t$. After a while, exposed hosts progress into the infectious but asymptomatic class, $A_t$. Then, they either enter the symptomatic class, $I_t$, or recover from the disease without showing symptoms. Eventually, symptomatic hosts either hopefully survive and recover or they die from the disease.

The total population is conserved over time, $S_t+E_t+A_t+I_t+R_t+D_t=1$, $t\ge 0$, and all processes but the infection one, are assumed to be based on geometric distributions, i.e., $\mathbb{P}(X=k)=p{(1-p)}^{k-1}$, $\mathbb{E}\left[X\right]=\frac{1}{p}$, $\mathrm{Var}\left(X\right)=\frac{1-p}{p^2}$, for some generic probability p. We recall that the geometric distribution describes the waiting time between successes in a Binomial process, i.e., when counting the number of successes in independent tries at something. Specifically, the model has five fixed probabilities: $0<\alpha ,\delta ,\gamma ,p,q<1$ corresponding to the mean latent period $\frac{1}{\alpha }$, the mean infectious period $\frac{1}{\delta }$ and $\frac{1}{\gamma }$ for the asymptomatic and symptomatic phases of infection, respectively, the probability of developing symptoms p, and the proportion of symptomatic hosts that result in death q.

Regarding to the non-linear infection process, we assume a force of infection, that is, the probability per unit of time of susceptible hosts being exposed to the pathogen, as the probability

$${\epsilon }_t=1-e^{-({\beta }_1{A}_t+{\beta }_2{I}_t)/(1-D_t)}$$

(1)

where ${\beta }_1,{\beta }_2>0$ are transmission rates in the asymptomatic and symptomatic phases of infection, respectively, and so depending on the number of infectious individuals over the alive population—see [8]. Actually, the expression in (1) corresponds to a Poisson probability distribution for the number of contacts per day between pathogen and susceptible host—see Chapter 9 in [6] for discrete-time models with a general force of infection.

According to the assumptions above, the discrete-time model equations for each epidemic outbreak reads as:

$$\left\{\begin{array}{c}{S}_{t+1}=(1-{\epsilon }_t)S_t\hfill \\ E_{t+1}={\epsilon }_t{S}_t+(1-\alpha )E_t\hfill \\ A_{t+1}=\alpha E_t+(1-\delta )A_t\hfill \\ I_{t+1}=p\delta A_t+(1-\gamma )I_t\hfill \\ R_{t+1}=(1-p)\delta A_t+(1-q)\gamma I_t+R_t\hfill \\ D_{t+1}=q\gamma I_t+D_t\hfill \end{array},t\ge 0,\right.$$

(2)

with suitable initial conditions at $t=0$.

For a better understanding of the epidemic model, let us recall some well-known definitions in regard to system (1) and (2):

• Incidence is the number of new cases per day in the interval $(t,t+1]$: ${\epsilon }_t{S}_{t}N=(S_t-S_{t+1})N$, where N would be the population size.
• Prevalence is the fraction of the population with infection or disease at a given time point: $E_t+A_t+I_t$.
• Force of infection ${\epsilon }_t$ is the probability per unit of time of a susceptible host becoming infected (here specifically, being exposed to the virus but not infectious yet).
• Infection transmission rates [1/day]: ${\beta }_1$, ${\beta }_2$ are contact rate × infectiveness probability. If $A_t+I_t\ll 1$, $D_t<1$, then we have ${\epsilon }_t=1-e^{-({\beta }_1{A}_t+{\beta }_2{I}_t)/(1-D_t)}\simeq {\beta }_1{A}_t+{\beta }_2{I}_t$.
• Mean incubation period, i.e., average time to symptom onset, is $\frac{1}{\alpha }+\frac{1}{\delta }$ days.
• Mean infectious period as either asymptomatic or symptomatic host is $\frac{1}{\delta }+\frac{1}{\gamma }$ days.
• Total number of cases until time t: ${\displaystyle \sum _{j=0}^{\infty }}{\epsilon }_{t-j}{S}_{t-j}N=(S_{-\infty }-S_{t+1})N$, where N would be the population size.
• Final size of the epidemics: fraction of the initial susceptible population that eventually becomes infected during the outbreak $(1-S_{\infty })$—see Section 6.
• Virulence, as a general concept [19,20], is the decrease in host fitness due to the infection. A measure of virulence is typically taken as the disease-induced mortality $q\gamma$, the severity of the pathogen.
• Case fatality ratio is the proportion of symptomatic cases that result in death q.
• Serial interval is the time between successive cases in a chain of transmission—see Section 4.

It is straightforward to check that system (1) and (2) has the disease-free steady state:

$$(S^{*},0,0,0,R^{*},D^{*}),\mathrm{with}{S}^{*}+R^{*}+D^{*}=1,$$

(3)

and as expected, no endemic equilibrium is possible because we have assumed no demographic turnover.

System (2) is inspired by the continuously structured model introduced in [24] where time since infection for the asymptomatic hosts is considered. However, the system of equations forms a discrete model for its own instead of being a simple discretization of a continuous model, since the force of infection ${ϵ}_t$ contains an exponential term instead of a truncated approximation—see the discussion in [8] and Chapter 9 in [6].

3. Non-Linear Renewal Equation

In order to reduce the model to a single non-linear renewal equation, we need to extend the initial conditions at $t=0$ to discrete histories in $(-\infty ,0]$, i.e., the state variables in the past, such that ${lim}_{j\to \infty }{S}_{-j}=1$, and ${lim}_{j\to \infty }{E}_{-j}={lim}_{j\to \infty }{A}_{-j}={lim}_{j\to \infty }{I}_{-j}={lim}_{j\to \infty }{R}_{-j}={lim}_{j\to \infty }{D}_{-j}=0$.

Firstly, we can reduce system (1) and (2) to 4 state variables since the fraction of removed and deceased hosts are given by

$$R_t=1-(S_t+E_t+A_t+I_t+D_t)\mathrm{and}{D}_t=q\gamma {\displaystyle \sum _{j=1}^{\infty }}{I}_{t-j},$$

(4)

respectively. Then, using the model equations recursively we get to:

$$\left\{\begin{array}{c}{S}_t={\displaystyle \prod _{j=1}^{\infty }}(1-{\epsilon }_{t-j})=exp\left(-\sum _{j=1}^{\infty }\frac{{\beta }_1{A}_{t-j}+{\beta }_2{I}_{t-j}}{1-D_{t-j}}\right)\hfill \\ E_t={\displaystyle \sum _{j=1}^{\infty }}{(1-\alpha )}^{j-1}{\epsilon }_{t-j}{S}_{t-j}\hfill \\ A_t=\alpha {\displaystyle \sum _{j=1}^{\infty }}{(1-\delta )}^{j-1}{E}_{t-j}\hfill \\ I_t=p\delta {\displaystyle \sum _{j=1}^{\infty }}{(1-\gamma )}^{j-1}{A}_{t-j}\hfill \end{array}\right.,t\ge 1,$$

(5)

Finally, plugging the equations one into another (formulas extended to any time t: past, present or future), we end up with the following scalar non-linear discrete renewal equation for the fraction of asymptomatic hosts $A_t$:

$$A_t=\alpha {\displaystyle \sum _{j=1}^{\infty }}{(1-\delta )}^{j-1}{\displaystyle \sum _{k=1}^{\infty }}{(1-\alpha )}^{k-1}{\epsilon }_{t-j-k}{\displaystyle \prod _{n=1}^{\infty }}(1-{\epsilon }_{t-j-k-n})$$

(6)

with ${\epsilon }_t=1-exp\left(-({\beta }_1{A}_t+{\beta }_{2}p\delta {\sum }_{j=1}^{\infty }{(1-\gamma )}^{j-1}{A}_{t-j})/(1-D_t)\right)$ and $D_t=pq\delta \gamma {\sum }_{k=1}^{\infty }{(1-\gamma )}^{k-1}{\sum }_{j=1}^{\infty }{A}_{t-j-k}$. Even though, non-linear renewal Equation (6) seems to be awkward, it has a clear and easy probabilistic interpretation. Indeed, rearranging the terms in (6), the fraction of asymptomatic hosts at time t is given by adding up the product of the following five terms:

$$\begin{array}{cc}\hfill A_t=\sum _{j,k\ge 1}\left[\prod _{n\ge 1}(1-{\epsilon }_{t-j-k-n})\right]·{\epsilon }_{t-j-k}·\alpha {(1-\alpha )}^{k-1}·\delta {(1-\delta )}^{j-1}·\frac{1}{\delta }=\sum _{j,k\ge 1}& \\ \hfill \mathrm{probability}\mathrm{for}\mathrm{a}\mathrm{host}\mathrm{of}\mathrm{being}\mathrm{susceptible}\mathrm{until}\mathrm{time}t-j-k& \times \hfill \\ \hfill \mathrm{probability}\mathrm{per}\text{time-unit}\mathrm{of}\mathrm{becoming}\mathrm{infected}\mathrm{at}t-j-k& \times \hfill \\ \hfill \mathrm{probability}\mathrm{latent}\mathrm{period}\mathrm{is}k\mathrm{days}& \times \hfill \\ \hfill \mathrm{probability}\mathrm{infectious}\mathrm{asymptomatic}\mathrm{period}\mathrm{is}j\mathrm{days}& \times \hfill \\ \hfill \mathrm{mean}\mathrm{infectious}\mathrm{asymptomatic}\mathrm{period}.& \end{array}$$

At this point we can compute the progression of the infection over time, either from the six evolution Equations (1) and (2) or from the single non-linear renewal Equation (6). On the one hand, we can compute the fraction of hosts in each disease stage from the model equations, taking suitable initial conditions at $t=0$ like

$$(S_0,E_0\simeq 0,A_0\simeq 0,I_0\simeq 0,R_0=0,D_0=0),\mathrm{with}{S}_0+E_0+A_0+I_0=1.$$

On the other hand, we can also compute the fraction of asymptomatic hosts $A_t$, $t\ge 1$, from the renewal equation (6) assuming that its initial history

$$\left\{A_{-j}\right|j\ge 0,\mathrm{with}\underset{j\to \infty }{lim}{A}_{-j}=0\}$$

is known. Then, we can compute the fraction of symptomatic hosts at time $t\ge 1$ as $I_t=p\delta {\displaystyle \sum _{j=1}^{\infty }}{(1-\gamma )}^{j-1}{A}_{t-j}$, with compatible initial condition $I_0=p\delta {\displaystyle \sum _{j=1}^{\infty }}{(1-\gamma )}^{j-1}{A}_{-j}$. Similarly, we can also compute the other state variables if needed.

4. Basic Reproduction Number

Taking advantage of the reduced systems (5) and (6), in this section we are going to compute the basic reproduction number ${\mathcal{R}}_0$ which describes the initial phase of the infection, [18,22,23,25,26].

We can compute the basic reproduction number for the present model once we have decided what is understood as an infection event. Here we have two natural points of view:

1.

the infection event is meant as the exposition to the pathogen of a susceptible host becoming an asymptomatic individual.

2.

the infection event is meant as the onset of symptoms for a host who has been exposed to the pathogen in the past.

Taking the first viewpoint, it is enough to linearize the renewal Equation (6), instead of linearizing (1) and (2) around the disease-free equilibrium.

First of all, around the disease-free steady state (3), the force of infection becomes ${\epsilon }_t\simeq {\beta }_1{A}_t+{\beta }_{2}p\delta {\displaystyle \sum _{j=1}^{\infty }}{(1-\gamma )}^{j-1}{A}_{t-j}$ and from non-linear renewal Equation (6), computed in the previous section, we get to the following linear discrete renewal equation:

$$A_t={\displaystyle \sum _{j=1}^{\infty }}\delta {(1-\delta )}^{j-1}{\displaystyle \sum _{k=1}^{\infty }}\alpha {(1-\alpha )}^{k-1}\left(\frac{{\beta }_1}{\delta }{A}_{t-j-k}+\frac{{\beta }_{2}p}{\gamma }{\displaystyle \sum _{n=1}^{\infty }}\gamma {(1-\gamma )}^{n-1}{A}_{t-j-k-n}\right)$$

(7)

In this equation, the three geometric distributions of the model appear clearly, with probabilities $\delta ,\alpha$ and $\gamma$, respectively. In other words, $\alpha {(1-\alpha )}^{k-1}·\delta {(1-\delta )}^{j-1}$ is the probability that the asymptomatic serial interval is $k+j$ days and $\alpha {(1-\alpha )}^{k-1}·\delta {(1-\delta )}^{j-1}·\gamma {(1-\gamma )}^{n-1}$ is the probability that the pre-symptomatic serial interval is $k+j+n$ days and so, (7) gives the asymptomatic hosts at time t computed from the asymptomatic ones in the past $s<t$. As expected, the average time between successive cases, with or without showing symptoms are $\frac{1}{\alpha }+\frac{1}{\delta }+\frac{1}{\gamma }$ and $\frac{1}{\alpha }+\frac{1}{\delta }$ days, respectively. Taking the standard approach of the next-generation matrix/operator—see [22,23,24,27,28], from the 1-dimensional renewal Equation (7) we readily get the basic reproduction number as:

$${\mathcal{R}}_{0,a}=\frac{{\beta }_1}{\delta }+\frac{{\beta }_{2}p}{\gamma }.$$

(8)

It is interpreted as the expected secondary asymptomatic cases produced by an asymptomatic primary case, since it is computed from the renewal equation for the asymptomatic hosts. However, this expression could have been derived directly from the model ingredients since, from a asymptomatic primary case, ${\mathcal{R}}_{0,a}=\frac{{\beta }_1}{\delta }+\frac{{\beta }_{2}p}{\gamma }$ traces the number of infections during the asymptomatic phase $\frac{{\beta }_1}{\delta }$ plus, provided that the host develops symptoms ($p>0$), the number of infections during the symptomatic phase $\frac{{\beta }_2}{\gamma }p$.

Furthermore, we can compute the effective reproduction number, that is, the time-dependent analog to the basic reproduction number (8), see [6], as

$${\mathcal{R}}_{t,a}=\frac{1-e^{-{\beta }_1{A}_t/(1-D_t)}}{A_t}·\frac{S_t}{\delta }+p\frac{1-e^{-{\beta }_2{I}_t/(1-D_t)}}{I_t}·\frac{S_t}{\gamma }.$$

(9)

It is an indicator of the transmission potential at each day of the epidemic outbreak and is analogously computed from the model ingredients as the expected number of new cases produced by an asymptomatic case, provided that the situation at the t-th day is unchanged.

See Figure 2 for an illustration of the progression over time of the three types of infected hosts (exposed → asymptomatic → symptomatic), computed from (1) and (2) with suitable initial conditions, as well as the effective reproduction number ${\mathcal{R}}_{t,a}$ that starts at $t=0$ as the basic reproduction number.

As we have already said at the beginning of the section, different reproduction numbers can be defined depending on what is understood as the infection event of the disease—see the discussion in [22]. Before tackling the computation of the basic reproduction number from the second viewpoint stated above, let us consider the following 2-dimensional linear discrete renewal equation:

$$\left\{\begin{array}{c}{I}_t=p\delta {\displaystyle \sum _{j=1}^{\infty }}{(1-\gamma )}^{j-1}{A}_{t-j}\hfill \\ A_t={\displaystyle \sum _{j=1}^{\infty }}{(1-\delta )}^{j-1}{\displaystyle \sum _{k=1}^{\infty }}\alpha {(1-\alpha )}^{k-1}\left({\beta }_1{A}_{t-j-k}+{\beta }_2{I}_{t-j-k}\right)\hfill \end{array}\right.$$

(10)

for the symptomatic and asymptomatic hosts $(I_t,A_t)$ (the ensemble of infectious hosts). It comes from last equation in (5), and (7). From system (10) we can get a joint basic reproduction number as the spectral radius of the 2-dim. matrix: $\left(\begin{array}{cc}0& \frac{p\delta }{\gamma }\\ \frac{{\beta }_2}{\delta }& \frac{{\beta }_1}{\delta }\end{array}\right)$ which turns out to be ${\tilde{\mathcal{R}}}_0=\frac{{\beta }_1}{2\delta }+\sqrt{{\left(\frac{{\beta }_1}{2\delta }\right)}^2+\frac{{\beta }_{2}p}{\gamma }}$ but it lacks of a clear epidemiological interpretation. However, we can use (10) to obtain a renewal equation for the symptomatic hosts. We proceed as follows. On the space of sequences (initial histories), we can define the linear operator ${\left(\mathcal{K}\varphi \right)}_t:={\sum }_{j=1}^{\infty }{(1-\delta )}^{j-1}{\sum }_{k=1}^{\infty }\alpha {(1-\alpha )}^{k-1}{\varphi }_{t-j-k}$. Then, the second equation in (10) becomes $A_t={\beta }_1{\left(\mathcal{K}A\right)}_t+{\beta }_2{\left(\mathcal{K}I\right)}_t$. If $\frac{{\beta }_1}{\delta }<1$, we can isolate $A_t={\beta }_2{\left({(Id-{\beta }_1\mathcal{K})}^{-1}\mathcal{K}I\right)}_t$ and (10) reduces to the following single renewal equation for $I_t$:

$$I_t={\beta }_{2}p\delta {\displaystyle \sum _{j=1}^{\infty }}{(1-\gamma )}^{j-1}{\left({(Id-{\beta }_1\mathcal{K})}^{-1}\mathcal{K}I\right)}_{t-j}.$$

Although we do not have the explicit expression of the inverse of the operator $\mathcal{K}$, again, from this 1-dimensional linear discrete renewal equation we readily get the basic reproduction number as:

$${\mathcal{R}}_{0,s}=\frac{{\beta }_{2}p}{\gamma }\frac{1}{1-{\beta }_1/\delta }\mathrm{whenever}\frac{{\beta }_1}{\delta }<1$$

(11)

where we have used that constant sequences are the eigenvectors of operators $\mathcal{K}$ and ${\mathcal{K}}^{-1}$. The expression above is interpreted as the expected number of symptomatic individuals that a symptomatic individual will produce. Taking into account that ${\mathcal{R}}_{0,s}=\frac{{\beta }_{2}p}{\gamma }\frac{1}{1-{\beta }_1/\delta }=\frac{{\beta }_2}{\gamma }p{\displaystyle \sum _{n=1}^{\infty }}{\left(\frac{{\beta }_1}{\delta }\right)}^{n-1}$, its direct derivation form the model ingredients would be as follows: from a symptomatic primary case, ${\mathcal{R}}_{0,s}$ is the sum of all pre-symptomatic cases generated from that primary case in their symptomatic phase: $\frac{{\beta }_2}{\gamma }\left(p+p\frac{{\beta }_1}{\delta }+p{\left(\frac{{\beta }_1}{\delta }\right)}^2+\cdots \right)$, provided that $\frac{{\beta }_1}{\delta }<1$. On the contrary, if $\frac{{\beta }_1}{\delta }\ge 1$ then the quantity ${\mathcal{R}}_{0,s}$ would not be bounded. We recall that p is the probability of developing symptoms.

As expected, the three expressions of the basic reproduction number: ${\mathcal{R}}_{0,a}$ in (8) from the asymptomatic viewpoint, ${\mathcal{R}}_{0,s}$ in (11) from the symptomatic viewpoint and ${\tilde{\mathcal{R}}}_0$ above from the joint viewpoint, are such that

$$\mathrm{sign}({\mathcal{R}}_{0,a}-1)=\mathrm{sign}({\mathcal{R}}_{0,s}-1)=\mathrm{sign}({\tilde{\mathcal{R}}}_0-1),$$

and they are related via a function of $\frac{{\beta }_1}{\delta }$ and $\frac{{\beta }_{2}p}{\gamma }$ and are independent of probabilities $\alpha ,q$. Let us remark that the expression in (8) is the one that the reader would expect intuitively, as it is for many reasons actually.

4.1. Initial Growth Rate

For the sake of completeness, let us determine the growth rate of the initial phase of the outbreak and check its expected relationship with the basic reproduction number of the infection. Let the initial growth rate of the disease be denoted by $\rho >0$. Looking for solutions of (7) as geometric progressions $A_t={\rho }_0·{\rho }^t$, we get to the following characteristic equation for $\rho >{\rho }_m:=max\{1-\delta ,1-\alpha ,1-\gamma \}$:

$$1=\frac{\alpha }{(\rho +\delta -1)(\rho +\alpha -1)}\left({\beta }_1+\frac{{\beta }_{2}p\delta }{\rho +\gamma -1}\right),$$

where $0<\alpha ,\delta ,\gamma ,p<1$ are four probabilities of the model and ${\beta }_1,{\beta }_2>0$ are transmission rates. The right hand side of this equation, as a function of $\rho >{\rho }_m$ is a continuous strictly decreasing function denoted by $F\left(\rho \right)$ such that $F\left(1\right)={\mathcal{R}}_{0,a}$, the basic reproduction number given in (8), ${lim}_{\rho \to {\rho }_m^{+}}F\left(\rho \right)=+\infty$ and ${lim}_{\rho \to \infty }F\left(\rho \right)=0$. Therefore, there is a unique real solution ${\rho }^{*}>{\rho }_m$ and, as expected, $\mathrm{sign}({\rho }^{*}-1)=\mathrm{sign}({\mathcal{R}}_{0,a}-1)$. Unfortunately, we cannot get to an explicit expression for the initial growth rate ${\rho }^{*}$.

So far all the computations are independent of the case fatality ratio q, except renewal equation (6) which is equivalent to the original Markov chain model.

5. Evolution of Infection Transmission

Once we have determined the occurrence of an epidemic outbreak, e.g., when ${\mathcal{R}}_{0,a}>1$, in this section we study the evolution of infection transmission before and after symptom onset. For similar studies on evolutionary epidemiology—see, e.g., [11,12,13,14,17,19,20,21].

We will focus on the long-term evolutionary dynamics of the virulence of the pathogen assuming some correlation between parameters of the model and taking for granted that evolution favors traits maximizing the basic reproduction number.

First of all, let us define the following weighted mean transmission rate:

$$\overline{\beta }=\frac{{\beta }_1+{\beta }_{2}p}{1+p}$$

(12)

which is an average between pre-symptomatic and post-symptomatic rate of infections ${\beta }_1$ and ${\beta }_2$, respectively, depending on the probability p of developing symptoms. If p would be zero then this mean would correspond to purely asymptomatic transmission.

Secondly, we will assume that virulence $q\gamma$, i.e., disease-induced mortality, is positively correlated with the mean transmission rate $\overline{\beta }$ during both infectious stages. Specifically, let us assume the following virulence-transmission tradeoff:

• Provided that hosts can develop symptoms $p>0$ and die from the disease $q>0$, virulence is positively correlated with the mean transmission rate (12) in the form:

  $$q\gamma =p·c{\overline{\beta }}^2\le 1$$

  (13)

  with $c>0$ a proportionality constant. Notice that if there was no possibility of developing symptoms then virulence would be null (no disease-induced mortality).

In this scenario, when mean transmission rate $\overline{\beta }$ increases, also virulence $q\gamma$ does according to (13) but transmission time (expected duration of both infectious stages) decreases since $T=\frac{1}{\delta }+\frac{1}{\gamma }=\frac{1}{\delta }+\frac{q}{pc{\overline{\beta }}^2}$. On the other hand, the fitness of the pathogen is typically measured by the basic reproduction number whose three equivalent expressions have been computed in Section 4 as:

$${\mathcal{R}}_{0,a}=\frac{{\beta }_1}{\delta }+\frac{{\beta }_{2}p}{\gamma },{\mathcal{R}}_{0,s}=\frac{{\beta }_{2}p}{\gamma }\frac{1}{1-{\beta }_1/\delta },{\tilde{\mathcal{R}}}_0=\frac{{\beta }_1}{2\delta }+\sqrt{{\left(\frac{{\beta }_1}{2\delta }\right)}^2+\frac{{\beta }_{2}p}{\gamma }}.$$

(14)

Without loss of generality, we can use the first expression (asymptomatic viewpoint) to optimize with respect to the transmission rate in the symptomatic phase. Indeed, putting (13) and (12) into ${\mathcal{R}}_{0,a}$ we get to the following function:

$${\mathcal{R}}_{0,a}\left({\beta }_2\right)=\frac{{\beta }_1}{\delta }+\frac{{\beta }_{2}q}{c{\overline{\beta }}^2}=\frac{{\beta }_1}{\delta }+\frac{{\beta }_{2}q}{c}{\left(\frac{1+p}{{\beta }_1+{\beta }_{2}p}\right)}^2$$

which has a global maximum at ${\beta }_2^{*}=\frac{{\beta }_1}{p}$, $p>0$, meaning that transmission rate is always higher after symptom onset ${\beta }_2^{*}>{\beta }_2^{*}p={\beta }_1$—see Figure 3 for an illustration with plausible values of model parameters. Analogous functions are obtained for ${\mathcal{R}}_{0,s}\left({\beta }_2\right)$ and ${\tilde{\mathcal{R}}}_0\left({\beta }_2\right)$.

Moreover, this optimal transmission rate can be translated into optimal virulence. Indeed, using the relationships in (12) and (13), we get to the optimal mean transmission ${\overline{\beta }}^{*}=\frac{2{\beta }_1}{1+p}$ and, most importantly, the optimal virulence $q{\gamma }^{*}=p·c{\left({\overline{\beta }}^{*}\right)}^2=p·c{\left(\frac{2{\beta }_1}{1+p}\right)}^2$. Therefore, under the virulence-transmission tradeoff of this section, there is an intermediate virulence for a pathogen to maximize the basic reproduction number—see Figure 4 for an illustration showing the optimal virulence (quite small level ~$8\%$) maximizing the basic reproduction number.

At optimal virulence $q{\gamma }^{*}$, since ${\mathcal{R}}_{0,a}^{*}={\beta }_1\left(\frac{1}{\delta }+\frac{1}{{\gamma }^{*}}\right)$ we have that most of the infections take place prior to symptom onset if $\delta <{\gamma }^{*}$, in other words, if the asymptomatic stage is longer than the symptomatic one—see Figure 5 for a summary of the outcomes of this section. See [15] for empirical data on COVID-19 suggesting that pre-symptomatic transmission may occur in a large proportion of transmission events.

6. Final Size of Symptomatic Hosts

In this section, we are interested in computing the final size of the symptomatic host population, i.e., the fraction of hosts who showed up symptoms at some point in the past. To do that, we have to compute the final size of the epidemics $1-S_{\infty }$ firstly. Let us remark that the computation below is rather technical.

The idea is to write the state variables depending on $S_{t-j}$, so in terms of the fraction of susceptible hosts in the past. First of all we have that ${\epsilon }_t=1-\frac{S_{t+1}}{S_t}$ from the first equation in (2). Then from (5),

$$E_t={\displaystyle \sum _{j=1}^{\infty }}{(1-\alpha )}^{j-1}(S_{t-j}-S_{t-j+1})\mathrm{and}{A}_t=\alpha {\displaystyle \sum _{j=1}^{\infty }}{(1-\delta )}^{j-1}{\displaystyle \sum _{k=1}^{\infty }}{(1-\alpha )}^{k-1}(S_{t-j-k}-S_{t-j-k+1}).$$

Now summing up the terms we get ${\displaystyle \sum _{n=1}^{\infty }}{A}_{t-n}=\alpha {\displaystyle \sum _{j=1}^{\infty }}{(1-\delta )}^{j-1}{\displaystyle \sum _{k=1}^{\infty }}{(1-\alpha )}^{k-1}(1-S_{t-j-k})$ and ${\displaystyle \sum _{n=1}^{\infty }}{I}_{t-n}=p\delta \alpha {\displaystyle \sum _{m=1}^{\infty }}{(1-\gamma )}^{m-1}{\displaystyle \sum _{j=1}^{\infty }}{(1-\delta )}^{j-1}{\displaystyle \sum _{k=1}^{\infty }}{(1-\alpha )}^{k-1}(1-S_{t-m-j-k})$. Finally, taking the limit as time goes to infinity, we have that ${\displaystyle \underset{t\to \infty }{lim}{D}_t=\underset{t\to \infty }{lim}q\gamma {\sum }_{n=1}^{\infty }{I}_{t-n}=pq(1-S_{\infty })}$ and we get to the “final” state:

$$S_{\infty }>0,E_{\infty }=0,A_{\infty }=0,I_{\infty }=0,R_{\infty }=(1-pq)(1-S_{\infty }),D_{\infty }=pq(1-S_{\infty }).$$

To end up, $S_t$ from the first equation in (5) is readily bounded in the interval:

$$exp\left(-{\displaystyle \sum _{n=1}^{\infty }}\frac{{\beta }_1{A}_{t-n}+{\beta }_2{I}_{t-n}}{1-D_{\infty }}\right)\le S_t\le exp\left(-{\displaystyle \sum _{n=1}^{\infty }}{\beta }_1{A}_{t-n}+{\beta }_2{I}_{t-n}\right)$$

and using the computations above, it turns out that ${\displaystyle \underset{t\to \infty }{lim}}{\displaystyle \sum _{n=1}^{\infty }}{\beta }_1{A}_{t-n}+{\beta }_2{I}_{t-n}=\left(\frac{{\beta }_1}{\delta }+\frac{{\beta }_{2}p}{\gamma }\right)(1-S_{\infty })$. Finally, we get a bounded interval for $S_{\infty }$ solving 2 non-linear equations:

$$e^{-{\mathcal{R}}_{0,a}\frac{1-S_{\infty }}{1-pq(1-S_{\infty })}}\le S_{\infty }\le e^{-{\mathcal{R}}_{0,a}(1-S_{\infty })}.$$

If $pq\ll 1$ we recover the classical equation $S_{\infty }=e^{-{\mathcal{R}}_{0,a}(1-S_{\infty })}$ for the final size of the epidemics—see [8]. Note again that ${\mathcal{R}}_{0,a}$ has a central role over the other expressions for the basic reproduction number.

Once we have determined the final size of the epidemics, either as the midpoint of the interval above for instance or computed by numerical simulation, we can compute the final size of the symptomatic hosts as $\frac{\gamma }{(1-p)\delta +\gamma }(1-S_{\infty })$. The latter takes into account the fraction of infected hosts reaching the state Removed or Deceased that corresponds to the symptomatic hosts—see the arrows in the diagram of Figure 1.

Taking parameter values of Figure 5 and using midpoint for the computation of $S_{\infty }$, we get that the final size of the epidemics is $(1-S_{\infty })=35.3\%$ and the final size of the symptomatic hosts is $\frac{\gamma }{(1-p)\delta +\gamma }(1-S_{\infty })=84.6\%(1-S_{\infty })=29.9\%$.

7. Generalization

If the goal is to fit the discrete-time model (2) to epidemiological data on a specific disease, some extensions are needed. Among others, for instance, we may wonder to include some loss of immunity along the same epidemic outbreak (i.e., allowing a fraction of recovered hosts to come back to the susceptible stage with reinfection probability $\theta$) and to extend the time at the infected stages $E\to A\to I$, from geometric distributions (discrete analog of exponential distributions) to negative binomial distributions (discrete analog of Gamma distributions). So, the Markov chain model (2) can be extended as the following enhanced model:

$$\left\{\begin{array}{c}{S}_{t+1}=(1-{\epsilon }_t)S_t+\theta R_t\hfill \\ E_{t+1}^1={\epsilon }_t{S}_t+(1-\alpha )E_t^1,E_{t+1}^i=\alpha E_t^{i-1}+(1-\alpha )E_t^i\hfill \\ A_{t+1}^1=\alpha E_t^n+(1-\delta )A_t^1,A_{t+1}^i=\delta A_t^{i-1}+(1-\delta )A_t^i\hfill \\ I_{t+1}^1=p\delta A_t^n+(1-\gamma )I_t^1,I_{t+1}^i=\gamma I_t^{i-1}+(1-\gamma )I_t^i\hfill \\ R_{t+1}=(1-p)\delta A_t^n+(1-q)\gamma I_t^n+(1-\theta )R_t\hfill \\ D_{t+1}=q\gamma I_t^n+D_t\hfill \end{array}\right.i=2\cdots n$$

(15)

where ${\epsilon }_t=1-e^{-({\beta }_1{A}_t+{\beta }_2{I}_t)/(1-D_t)}$ with $A_t={\sum }_{i=1}^n{A}_t^i$ and $I_t={\sum }_{i=1}^n{I}_t^i$ being here the total fraction of asymptomatic and symptomatic hosts, respectively, that is, the ensemble of hosts in the n sub-stages for each type.

Here, the underneath probabilistic model for the linear transitions between infected stages $E\to A\to I$, is given by $\mathbb{P}(X=k)=\left(\genfrac{}{}{0pt}{}{k-1}{n-1}\right)p^n{(1-p)}^{k-n}$, $k\ge n$, $\mathbb{E}\left[X\right]=\frac{n}{p}$, $\mathrm{Var}\left(X\right)=n\frac{1-p}{p^2}$, for some generic probability p and some fix integer $n\ge 1$. As in Section 2, the total population is conserved over time and the model has six fixed probabilities: $0<\alpha ,\delta ,\gamma ,p,q,\theta <1$

For system (15), one can check that only the disease-free steady state is possible $E_{*}^i=A_{*}^i=I_{*}^i=0$, $i=1\cdots n$, $S^{*}+R^{*}=1$ and $D^{*}=0$, and the basic reproduction number from the asymptomatic viewpoint is analogous to (8) and reads as:

$${\mathcal{R}}_{0,a}=\frac{{\beta }_{1}n}{\delta }+\frac{{\beta }_{2}pn}{\gamma },$$

where we recall that n is the number of sub-stages considered. Again, this expression can be readily derived from the model ingredients taking into account that here the expected asymptomatic/symptomatic infectious period is $\frac{n}{\delta }$ and $\frac{n}{\gamma }$, respectively—see also Section 4.

Analogous equations to (6) and (7) can be derived for the enhanced model (15), either by the probabilistic interpretation (negative binomial distributions) or by the computation of the terms involved.

8. Discussion and Conclusions

Discrete-time epidemic models are simple yet powerful dynamical systems to describe the spread of an infectious disease. They are very easy to implement on a computer. For the present model (1) and (2) in particular, we have achieved a probabilistic interpretation of the model when reducing the whole system to a single non-linear renewal equation for the fraction of asymptomatic hosts. See the detailed interpretations of Equations (6) and (7), in particular, the determination of the probability distribution of serial intervals (time between successive cases, with or without showing symptoms). The importance of renewal equations in ecology and epidemiology has been already identified by many authors in the past. See the books by N. Bacaër et al. [29,30,31] and the references therein for a review on the history of models of population dynamics. In addition, we have focused on two important indicators, namely, the transmission potential (basic reproduction number) and the severity of the pathogen (virulence).

For a model with two different type of infectious hosts, the ones showing mild or no symptoms and the ones showing clear symptoms of the disease, it is interesting to wonder which viewpoint is most important. According to our findings, the asymptomatic viewpoint is the most important when reducing the model to a renewal equation in Section 3 or computing the final size of the epidemics in Section 6, or when deriving the basic reproduction number in Section 4 (and the related non-explicit initial growth rate $\rho$ of the infection in Section 4.1). The expression of ${\mathcal{R}}_{0,a}$ in (8) is what one would expect intuitively as the basic reproduction number, i.e., the sum of the transmission rate times the expected duration of the transmission in each infectious stage, taking into account the proportion of hosts reaching the symptomatic phase. Its counterpart from the symptomatic viewpoint is mathematically more involved, and most importantly, we have to convince the reader that the sought expression (11) corresponds to a properly-defined basic reproduction number actually. Regarding the infection transmission before and after the symptom onset, the best strategy for a pathogen is to silently spread as much as possible until the symptom onset and then enhance the transmission rate to compensate the fact that some of the hosts will not reach the symptomatic phase (i.e., they will recover without showing symptoms). The approach and the assumptions taken in Section 5 allows us to draw this conclusion—see Figure 5 summarizing the outcomes.