Abstract
We propose two variants of a stochastic epidemic model in which the disease is spread by mobile particles performing random walks on the complete graph. For the first model, we study the effect on the epidemic size of an immunization mechanism that depends on the activity of the disease. In the second model, the transmission agents can gain lives at random during their existences. We prove limit theorems for the final outcome of these processes. The epidemic model with mutations exhibits phase transition, meaning that if the mutation parameter is sufficiently large, then asymptotically all the individuals in the population are infected.
1. Introduction
One of the most important applications of stochastic processes is in epidemiology, more specifically, in studying the spread of epidemics and analyzing mechanisms to control them. The first studies of stochastic systems of epidemics go back to a pioneer paper of McKendrick [1] and the unpublished lecture notes of Reed and Frost, written around 1930. But only since the late 1940s, with Bartlett [2], that stochastic epidemic models began to be studied more extensively. Recently, the mathematical theory of epidemics has experienced a renewal impulse, due to the extension of the concept of a disease, to include computer malware, rumour spreading, and viral marketing. For well-crafted presentations of this theory, the reader should refer to the books by Anderson and May [3], Andersson and Britton [4], Bailey [5], Daley and Gani [6], and Draief and Massoulié [7].
In many models of the literature on epidemics, individuals reside at the vertices of a graph and the disease is transmitted through contacts along the edges, in such a way that the state of an individual depends on the state of other individuals in its vicinity. Typically, one considers the complete graph (representing a homogeneously mixing population) and the dissemination of the disease is described by a continuous-time Markov process whose transition probabilities are nonlinear in nature. The classical approach for studying such stochastic epidemic processes is to analyze the forward equation, which frequently results in differential equations that are very hard to tract analytically. In the last decades, the main advances have been made by using other tools, such as the analysis of the embedded Markov chain, coupling, martingale arguments, and diffusion approximations. In our paper, we study the final outcome of two epidemic processes described through interacting particle systems. Although our models have the standard assumption of homogeneous mixing, the possible interactions are dynamic: we suppose that the viruses are mobile particles which move along the vertices of the graph. In addition, the classical coupling techniques are difficult to apply, by virtue of the dependence between the progeny and the lifetime of the pathogenic agents. Thus, we illustrate how the machinery of density dependent population processes is useful for proving results about the final outcome of stochastic epidemic models on homogeneously mixing populations. This theory, described in Ethier and Kurtz [8, Chapter 11], is also exploited for this purpose in Ball and Britton [9, 10] and Kurtz et al. [11]. A detailed treatment of classical epidemic models on general networks can be found in Draief and Massoulié [7, Chapter 8].
Let us describe the dynamics of the interacting particle system formulated in Kurtz et al. [11], which models the propagation of a disease in a population and the spreading of a virus in a computer network. For , let be the complete graph on vertices, that is, the undirected graph with vertices where each pair of vertices is connected by an edge. At time zero, there is one particle at each vertex of ; one of them is active; the others are inactive. Active particles perform independent, continuous-time, rate 1, random walks on , and, whenever any active particle visits an inactive one, the latter becomes active and starts its own random walk. However each active particle dies at the instant it jumps on a vertex that has already been visited by the process.
To understand the epidemic interpretation, think that the active particles are viruses which move along the vertices of the graph (individuals or computers). During its life, at the time points of a homogeneous Poisson process with intensity 1, a virus moves to a randomly chosen vertex of the graph, different from the one at which it is located. When the virus jumps onto a susceptible individual (i.e., a vertex that has not already been visited by the process), this individual catches the disease and the virus duplicates. Once that happens, the individual is regarded as immunized; that is, it activates an antivirus which will kill any virus that tries to infect it in the future. Thus, we may think that each newly born virus starts with one life, and it dies whenever it jumps onto an immunized individual. Since the process eventually finishes (when all viruses die), the main question concerns the proportion of individuals who ultimately become infected (i.e., the proportion of visited vertices at the end of the process). Kurtz et al. [11] prove a Law of Large Numbers for this proportion, which states that, for large , approximately 80% of the individuals are infected with high probability. A Central Limit Theorem describing the magnitude of the random fluctuations around this limiting value is also established.
Our purpose is to generalize these limit theorems for two variants of this epidemic model in which we introduce mechanisms of vaccination and mutations of the viruses. In the first model, presented in Section 2, we incorporate a process of vaccination in the epidemic model that depends on the activity of the disease and analyze how this affects the asymptotic distribution of the final epidemic size. In particular, we are able to evaluate the immunization rate required to keep the epidemic at a reasonable level with high probability. In the second model, described in Section 3, we consider a mutation parameter , which determines the epochs when the viruses gain more lives. This model exhibits a phase transition in the parameter, which means that for asymptotically the infection is restricted to a fraction less than of the population, whereas for asymptotically all the individuals are infected. In Section 4, we present some numerical illustrations of the main results for both models. Section 5 is devoted to the proofs.
2. Epidemic Model with Immunization
Let us first review the basic definitions and main results stated for the epidemic model formulated by Kurtz et al. [11]. For each and , let and denote, respectively, the number of visited vertices and the number of active particles at time , where . Then, the evolution of the epidemic is described by a continuous-time Markov chain in which proceeds according to the following transition scheme: The two possible transitions and correspond, respectively, to the occurrence of a contagion (infection of a susceptible individual and creation of a new virus) and the occurrence of a removal (death of a virus). The transition rates have the usual meaning; for instance, the first one means that
Let be the time until the epidemic ceases. Notice that if all individuals are infected then the process evolves like a death chain in the second coordinate, so that is finite almost surely. Thus, the proportion is the final epidemic size (which is called the coverage of in Kurtz et al. [11]). To state the limit theorems for this proportion, let denote the principal branch of the Lambert function. More details on this function may be found in the Appendix. From the results proved in Kurtz et al. [11], we have that where Furthermore, where denotes convergence in distribution and is the Gaussian distribution with mean zero and variance given by This means that, for large , approximately 80% of the individuals are infected with high probability. In addition, the distribution of the final epidemic size is asymptotically normal with mean and variance .
Now we define the first variant of the epidemic model, in which we introduce an infection control mechanism, like an immunization treatment or vaccination. We consider the possibility of immunizing the individuals of the population (vertices of the graph), in order to reduce the epidemic outbreak, and study how this influences the final outcome of the process. We suppose that we cannot kill the active particles (viruses) directly and incorporate a process of vaccination in the epidemic model that is tied to the intensity of the epidemic as it evolves. We assume that, apart from the natural dynamics of the model, at rate proportional to the number of active particles, a vertex is chosen at random and is immunized if it has not been visited by the process (i.e., if it is still susceptible). Whenever an active particle jumps onto an immunized vertex, it dies. This corresponds to a door-to-door immunization process, which may occur at each time step of the embedded Markov chain. We denote by the immunization rate and let be the number of immunized individuals at time . Then, the process is a continuous-time Markov chain that starts at the point and makes three possible transitions during each time increment : These possible transitions correspond, respectively, to the occurrence of a contagion, a removal of a virus, and an immunization of a susceptible individual. The associated infinitesimal transition rates are as follows: Let denote the time until the epidemic process stops (as given in (3)), so the proportions of infected and immunized persons in the population after all viruses die are given, respectively, by
Our main results are a Law of Large Numbers and a Central Limit Theorem for these proportions.
Theorem 1. Define Then,
Theorem 2. One has that where is the bivariate normal distribution with mean zero and covariance matrix given by
The proofs of Theorems 1 and 2 are given in Section 5. The dependence of the limiting fractions of infected and immunized individuals on is illustrated in Figure 1(a). As one would expect, increasing the immunization rate pushes the value of towards zero. This means that the epidemic outbreak can be controlled by a suitable choice of the immunization rate. In Figure 1(b), the components of the covariance matrix are plotted as a function of .
(a)
(b)
We underline that the Central Limit Theorem allows us to evaluate the immunization rate required to keep the epidemic at a reasonable level with high probability. Indeed, given a large value of , and numbers and , we can compute the value of for which the probability that is approximately equal to . Table 1 exhibits the required values of for some arbitrarily chosen values of , , and .
3. Epidemic Model with Mutations
Now we define a version of the epidemic model with mutations, in which a virus is allowed a random number of jumps to already visited vertices before dying (i.e., it has a random number of lives). In biological terms, this corresponds to the ability of the pathogenic agents to develop resistance, by evolutionary adaptation. We consider the initial configuration with one particle per vertex of , only one of them being active. We assume that each newly activated particle starts with one life and assign to it two independent Poisson clocks, with rates 1 and . Poisson clocks assigned to different particles are independent of each other. The rate Poisson clock determines the mutation times: every time it rings, the active particle gains one life. As usual, the rate 1 Poisson clock determines the jump epochs of the particle. If the chosen vertex has already been visited, then the active particle loses one life (dying if it has one life); otherwise the inactive particle on this vertex is activated. We call the mutation rate.
Let be the number of visited vertices at time , and let be the number of active particles with lives at time (). We also define the number of active particles at time . Notice that the infinite-dimensional process makes transitions in continuous time according to the following transition rates: We underline an important feature of this epidemic model: it may continue forever, even if every individual is infected by the disease. We consider the process as finished either when there are no more active particles or when the whole population is infected; that is, we define its absorption time by Thus, is the final epidemic size.
Next we present the limit theorems for the final epidemic size. The Law of Large Numbers states that the epidemic model with mutations exhibits phase transition, in the sense that the limiting value of the final epidemic size is less than or equals accordingly as or . That is, if the pathogenic agents are able to gain resistance at a very large rate, then a global outbreak affecting all individuals is almost certain to occur.
Theorem 3. (i) For , one defines
Then,
(ii) If , then
Theorem 4. If , then where
The graphs of the limiting fraction of infected persons and the variance of the asymptotic normal distribution as functions of for are shown in Figure 2.
(a)
(b)
4. Numerical Illustrations
In order to illustrate the asymptotic theorems presented in Sections 2 and 3, we performed simulation studies of both models. Figure 3 is based on simulations of the epidemic model for a population consisting of individuals, with immunization parameter . Histograms (a) and (b) correspond to the rescaled values and obtained from the simulations, respectively. The empirical means of the final proportions of infected and immunized individuals were and , whereas the corresponding theoretical values are and , almost the same as the empirical means. The limiting covariance matrix has components , , and . The corresponding empirical values from the simulations were , , and , again close to their asymptotic counterparts. Overlaid on the histograms are the probability density functions of the limiting Gaussian distributions having zero mean and variances and . Figure 4 shows the ellipsoid quartiles of the rescaled data obtained from the simulations and the , , and contours of the corresponding asymptotic bivariate normal density. As can be seen, there is an excellent agreement between the empirical distributions and the theoretical limiting distributions.
(a)
(b)
We also ran simulations of the epidemic model with population size and mutation parameter . For this choice of , the theoretical asymptotic values are and . The corresponding estimates obtained from the simulations were and . Figure 5 presents the histogram of the rescaled final epidemic size from the simulations, with the limiting normal probability density function superimposed.
5. Proofs
A fundamental approach to study the behaviour of stochastic epidemic processes is to analyze their large population deterministic limits and diffusion approximations. In our proofs, we achieve that by applying the theory of density dependent population processes, described in Ethier and Kurtz [8, Chapter 11]. In particular, we will rely on Theorem 11.4.1 therein, which enables limit theorems to be derived for the final outcome of the epidemic models. To use this result, however, we have to first change the pace at which the process evolves, without affecting where it is absorbed. The main steps follow a path similar to that presented in Kurtz et al. [11].
5.1. Proofs of Theorems 1 and 2
For , we define From (9), we conclude that the continuous-time Markov chain has the following table of possible transitions and associated rates, all relating to a time increment : Thus, the transition rates depend on the current state of the process only through the density process , so is a density dependent population process. Observe that we are dealing with the epidemic model on , but the denominator in (22) is equal to . We adopt this definition only to facilitate the application of the related theory.
Let denote the absorption time of the process . Since Theorems 1 and 2 are established once we prove the following lemmas.
Lemma 5. and in probability, where and are given by (11).
Lemma 6. as , where is given by (14).
To prove these results, the first step is to define a new version of the process , for which we can apply Theorem 11.4.1 of Ethier and Kurtz [8].
5.1.1. Time-Changed Process
Consider the discrete-time embedded Markov chain of the process , that is, the discrete-time Markovian structure obtained from by ignoring the waiting times between consecutive jumps. It is not difficult to see that the distribution of depends on only through this embedded Markov chain. Hence, by defining the waiting times between consecutive transitions suitably, we can construct a coupled process with the same transitions as the original process until absorption, though with time at a different pace. This construction is done in such a way that the process has associated transition rates given by where Furthermore, defining , we have that
Now we define the functions so the infinitesimal rates in (26) admit the form . Consequently, is a density dependent population process with possible transitions in the set .
5.1.2. Deterministic Limit of the Time-Changed Process
Towards proving Lemmas 5 and 6, we use Theorem 11.2.1 of Ethier and Kurtz [8] to conclude that the time-changed system converges almost surely as (on a suitable probability space). The corresponding drift function is given by hence we obtain the following system of ordinary differential equations for the large population limit deterministic system: To solve this system, we define , . From (31), it follows that satisfies the ordinary differential equation with initial condition , whence Now we define the function given by Then, by solving each ordinary differential equation in (31), we obtain that the solution of this system is given by , where It follows from Theorem 11.2.1 of Ethier and Kurtz [8] that, as , converges almost surely to , uniformly on any finite time interval. Furthermore, analogously to Lemma 3.6 in Kurtz et al. [11], it can be proved that converges almost surely to , uniformly on , and that the same assertion holds for and .
5.1.3. Proofs of Lemmas 5 and 6
Both results follow from Theorem 11.4.1 of Ethier and Kurtz [8] applied to the time-changed process (recall formula (28)). We adopt the notations used there, except for the Gaussian process defined on page 458, that we would rather denote by . Here , and We also observe that given in (11) is the unique root of the function in the interval (see the Appendix for more details). In addition, the function in (32) establishes a one-to-one correspondence between and . Consequently, defining we conclude from (11) and (34) that , , and . Thus, to establish Lemma 5, it is enough to prove that The latter is exactly the first main statement of the proof of Theorem 11.4.1 from Ethier and Kurtz [8], whose argument we recall. The facts that and imply that and for . Then, the almost sure convergence of to uniformly on any finite time interval yields (37).
With regard to Lemma 6, we obtain from Theorem 11.4.1 of Ethier and Kurtz [8] that converges in distribution as to The resulting asymptotic distribution is a mean zero bivariate normal distribution, whose covariance matrix can be computed with the aid of a mathematical software. First, the matrix of partial derivatives of the drift function and the matrix are obtained as Furthermore, the solution of the matrix equation is given by Thus, using that the covariance matrix of the Gaussian process is we get that Finally, we obtain formula (14) by using the well-known properties of variance and covariance and simplifying properly.
5.2. Proofs of Theorems 3 and 4
For , we define Recalling (15), we have that the process makes transitions in continuous time according to the following transition rates: Denoting the absorption time of this process by we obtain that Consequently, Theorems 3 and 4 are established once we prove the following results.
Lemma 7. If , then in probability, and as , where and are given, respectively, by (17) and (21).
Lemma 8. If , then in probability.
The main idea to prove Lemma 7 is to define a time-changed version of the process and then apply Theorem 11.4.1 of Ethier and Kurtz [8] to a density dependent Markov chain in obtained from . The proof of Lemma 8 is similar, except for the fact that we work with another time-changed reduced Markov chain. Since the arguments are analogous to those presented in Section 5.1 and in Kurtz et al. [11], we present only the main steps of the proofs.
5.2.1. Proof of Lemma 7
We can define a coupled process with the same transitions as the process until absorption, though with time at a different pace. In view of (46), the process is constructed in such a way that it has the same possible transitions as during each time increment, but the associated rates are given by where We observe that corresponds to in the new time-scale. Moreover, defining , we have that . In order to prove Lemma 7 using Theorem 11.4.1 of Ethier and Kurtz [8], we work with a reduced Markov chain. We define and let , . Notice that is a Markov process that summarizes the evolution of the infinite-dimensional process . During each time increment , the possible state transitions associated to are , , and . In fact, we have that is a density dependent Markov chain with corresponding transition functions given by Lemma 7 follows from Theorem 11.4.1 of Ethier and Kurtz [8] applied to . Here, we have that The corresponding system of ordinary differential equations is To solve this system, we define the function given by . The solution is , where As , we have that converges almost surely to , uniformly on bounded time intervals (on a suitable probability space).
Notice now that given in (17) is the unique root of the function in . Moreover, the function in (55) establishes a one-to-one correspondence between and , so , where . Since and we conclude that almost surely. Consequently,
Regarding the Central Limit Theorem, the associated Gaussian process is denoted here by , and we obtain from Theorem 11.4.1 of Ethier and Kurtz [8] that The principal points in the calculation of the variance of the asymptotic normal distribution are the following: Again, by using the well-known properties of variance and covariance and simplifying properly, we obtain formula (21).
5.2.2. Proof of Lemma 8
Notice that we can construct the epidemic process with mutation rate from the process with rate , in such a way that the final number of infected individuals in the latter process is no greater than in the former. Consequently, it is enough to prove that in probability only for the process with rate . Considering such a process, let be as defined previously, and construct from a time-changed process which is a density dependent Markov chain with the same transitions as , but with transition functions given by Of course, this process is defined up to .
The solution of the corresponding system of ordinary differential equations is given by Since converges almost surely to uniformly on , we conclude that almost surely. From this, it follows that almost surely.
Appendix
The Lambert Function
The Lambert function is the multivalued inverse of the function . That is, is defined to be the function satisfying If is real, then for there are two possible real values of (see Figure 6(a)). The branch that satisfies is called the principal branch of the function and denoted by ; the branch that satisfies is called the lower branch and denoted by . We refer to Corless et al. [12] for more details.
(a)
(b)
Now let and consider the function given by The graph of is presented in Figure 6(b). It is not difficult to see that , is increasing in and decreasing in , and . Therefore, has a unique root in the interval .
To express in terms of the function, define , thus . Notice that satisfies , which can be written as Consequently,
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The author is grateful to CNPq (Grant 303872/2012-8), FAPESP (Grant 2012/22673-4), and PRP-FAEPEX-UNICAMP (Grant 016/13) for financial support. He also thanks the two anonymous referees for valuable suggestions which allowed improving the paper.