key: cord-0266471-r0q73a5f
authors: Cai, Chao-Ran; wu, Zhi-Xi; Guan, Jian-Yue
title: Behavior of susceptible-vaccinated--infected--recovered epidemics with diversity in the infection rate of the individuals
date: 2013-12-03
journal: nan
DOI: 10.1103/physreve.88.062805
sha: 4fcab79d5b0cb66dcbb4c9a49b16a2a37da331fc
doc_id: 266471
cord_uid: r0q73a5f

We study a susceptible-vaccinated--infected--recovered (SVIR) epidemic-spreading model with diversity of infection rate of the individuals. By means of analytical arguments as well as extensive computer simulations, we demonstrate that the heterogeneity in infection rate can either impede or accelerate the epidemic spreading, which depends on the amount of vaccinated individuals introduced in the population as well as the contact pattern among the individuals. Remarkably, as long as the individuals with different capability of acquiring the disease interact with unequal frequency, there always exist a cross point for the fraction of vaccinated, below which the diversity of infection rate hinders the epidemic spreading and above which expedites it. The overall results are robust to the SVIR dynamics defined on different population models; the possible applications of the results are discussed.

Infectious diseases have always been the great enemy of human health. Historically, large outbreaks of epidemic usually posed a great threat to health and caused great loss for individuals. In some sense, the history of humans is a history of struggle with all kinds of diseases, from the Black Death in medieval Europe to the recently notorious severe acute respiratory syndrome [1] [2] [3] , avian influenza [4, 5] , swine influenza [6, 7] , etc.

So far, vaccination is the most effective approach to preventing transmission of vaccine-preventable diseases, such as seasonal influenza and influenza like epidemics, as well as reducing morbidity and mortality [8] . In a voluntary vaccination program, the individuals are subject not only to social factors such as religious belief and human rights, but also to various other conditions such as risk of infection, prevalence of disease, coverage, and cost of vaccination.

Recently, a great deal of effort has been devoted to the investigation of the interplay between vaccine coverage, disease prevalence, and the vaccinating behavior of individuals by integrating game theory into traditional epidemiological models [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] . For brief reviews of this research topic, we refer the reader to Refs. [20, 21] and reference therein. Bauch et al. used game theory to explain the relationship between group interest and self-interest in smallpox vaccination policy [8, 9] and found that voluntary vaccination was unlikely to reach the group-optimal level. Vardavas and co-workers investigated the effect of voluntary vaccination on the prevalence of influenza based on minority game theory and showed that severe epidemics could not be prevented unless vaccination programs offer incentives [10] . Zhang et al. studied the epidemic spreading with voluntary vaccination strategy on both Erdös-Rényi random graphs and Barabási-Albert scalefree networks [12] . They found that disease outbreak can be more effectively inhibited on scale-free networks rather than * wuzhx@lzu.edu.cn † guanjy@lzu.edu.cn on random networks, which is attributed to the fact that the hub nodes of scale-free networks are more inclined to getting vaccinated after balancing the pros and cons. More recently, Fu and co-workers proposed a game-theoretic model to study the dynamics of vaccination behavior on lattice and complex networks [13, 15] . They found that the population structure causes both advantages and problems for public health: It can promote voluntary vaccination to high levels required for herd immunity when the cost for vaccination is sufficiently small, whereas small increases in the cost beyond a certain threshold will cause vaccination to plummet, and infection to rise, more dramatically than in well-mixed populations. Another research line studying the effect of human behavior on the dynamics of epidemic spreading considers mainly the coevolution of node dynamics and network structure (the so-called adaptive networks [22] ), which can affect considerably the spreading of a disease [23] [24] [25] .

In most classical epidemiological models [26, 27] , the individuals in the population are assumed to be identical, e.g., all susceptible individuals acquire the disease with the same probability whenever in contact with an infected individual, and all infected individuals recover, or go back to being susceptible, with the same rate. Such consideration is, however, far from the actual situation. Generally, catching a disease could be caused by many complex factors and there might be great difference among the individuals in the contact rate [28] , the infection rate (or disease transmission rate) [29, 30] , the recovery rate, the cost when the individual is infected, and so forth. One example of such a scenario would be the case where the population is divided into a relatively wealthy class (e.g., representing urban residents), which is less susceptible to infectious disease being considered due to better living conditions and/or health care, and a class of relatively impoverished (e.g., representing rural residents), which is more susceptible to infection. An alternative view is to regard roughly the whole population as composed of two main groups, say, youths and adults, where the former is more resistant to disease than the latter, owning to their stronger physique and immune system.

In the present work, we relax the assumption of identical nature of the individuals and take into account their hetero-geneity in acquiring disease when in contact with infectious individuals. To do this, we divide the whole population into two groups, youths (hereafter group A) and adults (group B) for simplicity, with the same size, and assume that the individuals from group B are more likely to be infected than those from group A. For the sake of comparison, we presume that only the disease transmission rate for the individuals in the two groups are distinct and other parameters, including the recovery rate, the cost of infection, and the cost for vaccination, are identical. By doing so we hope to catch any possible effects on disease prevalence and vaccination coverage caused by the variability of susceptibility. Our results presented below show that the heterogeneity in infection rate has a significant influence on disease spreading and hence cannot be ignored in the forecast of epidemic size and vaccination coverage.

Our paper is organized as follows. In Sec. II we define our model and give detailed information for the numerical simulation method and the parametrizations. In Sec. III we present and analyze the main results of our model. We summarize and discuss the results in Sec. IV.

It is well known that the contact pattern among individuals dramatically impacts the spatiotemporal dynamics of epidemic spreading in a population [26, 27] . In order to examine the robustness of the results of our model, we consider two types of population models, namely, a simple metapopulation model and a spatially structured population model, as illustrated in Fig. 1 .

In the metapopulation model, the whole population is divided into two subpopulations with equal size, namely, group A and group B. Within each subpopulation, the individuals are assumed to be homogeneously mixed, that is, every individual has the same opportunity to be in contact with everyone else. Generally speaking, because of the diversity in social conditions or lifestyles, the individuals living in an urban area would be more likely to interact with those also living the same area and less likely to interact with those in the suburb. Therefore, we consider the distinct contact pattern among the individuals to study its impact. This is done by assuming that any pair of individuals from different (the same) groups have an interaction frequency (1-). Here is restricted to the interval [0,0.5]. In the spatially structured population model, we consider two kinds of occupation of the individuals on a square lattice to introduce the diversity of interaction pattern among them. To be more specific, in the first case, the youths and the adults are arranged in a random way such that they can interact with the same frequency, which is similar to the case of = 0.5 in the metapopulation case. In the second case, the individuals are regularly prearranged to gather together with the same type of individuals [see Fig. 1(c) ]. In this way, we are able to investigate how the mixing pattern affects the epidemic spreading in the population.

We implement our susceptible-vaccinated-infectedrecovered epidemic-spreading dynamics in the following way. The epidemic strain infects an initial number of individuals I 0 and then spreads in the population according to the classical susceptible-infected-recovered(SIR) epidemiological model, with per-day transmission rate r for each pair of susceptible-infected contact and recovery rate g for each infected individual getting immune to the disease. Whenever the vaccinated compartment is involved in the epidemiological model, a fraction f V of individuals are randomly chosen in the whole population in the initial stage to get vaccinated. For simplicity, here we assume that vaccination grants perfect immunity for the infectious disease. The epidemic continues until there are no more newly infected individuals. As such, those unvaccinated susceptible individuals would either be infected or successfully escape from infection at the end of each spreading season.

In realistic situations, to vaccinate or not to vaccinate is sometimes the business of the individuals. Thus, except for the above case where the fraction of vaccinated individuals is compulsively introduced, we also consider a voluntary vaccination program for preventing an influenzalike infectious disease, in which individuals need to decide whether or not to receive a vaccine each season based on their perceived risk of disease infection. Following previous studies [10, 11, 13, 15] , we model the vaccination dynamics as a two-stage game. At the first stage, each individual decides whether or not to get vaccinated, which will incur a cost C V , including the immediate monetary cost for vaccine and the potential risk of vaccine side effects. Individuals catching the epidemic will suffer from an infection cost C I , which may account for disease complications, expenses for treatment, etc. Those individuals who escape infection are free riders and pay for nothing. Without loss of generality, we set C I = 1 and let c = C V /C I describe the relative cost of vaccination, whose value is restricted in the region of [0,1] (otherwise, doing nothing would be better than getting vaccinated). The second stage is the same epidemic spreading processes as described before. After each spreading season, the individuals are allowed to rechoose their choice for vaccination based on a pairwise comparison rule (more details will be given below).

We carry out stochastic simulations for the above epidemiological (game-theoretic) processes in both population models, wherein each seasonal epidemiological process is implemented by using the well-known Gillespie algorithm [31, 32] . In particular, at any time t, we calculate each individual's transition rate λ i (t). The rate for any susceptible individual becoming infected is λ i (t) = r × k inf and k inf is the number of infected neighbors of the focal individual. The rate for any infected individual recovering is λ i (t) = g. The recovered individuals do not change state and the rate for them is therefore λ i (t) = 0. Summing up all of them, we yield the total transition rate ω(t) = Σ i λ i (t). With this value in hand, the time at which the next transition event occurs is t = t + ∆t, where ∆t is sampled from an exponential distribution with mean 1 ω(t) (if we generate a uniform random number u ∈ [0, 1), then the time interval is ∆t = − ln(1−u) ω(t) ). The individual whose state is chosen to change at time t is sampled with a probability proportional to λ i (t). That is, 

, then individual k is chosen to change state. This elementary step is repeated until there are no infected individuals left in the population.

We first examine our model in metapopulations. For convenience, the two groups A and B are denoted by the subscripts a and b, respectively. According to the above illustrated scenario, the time evolution of population states for group A can be expressed as the following deterministic ordinary differential equations:

As mentioned before, the parameter is the cross contact coefficient, which stands for the contact frequency between individuals from different groups. For the whole system that includes groups A and B, we have the following equations:

In the limit → 0, the basic reproduction number (whose value identifies the expected number of secondary infections produced by an infected individual during that individual's infectious period within the entire susceptible population) of groups A and B can be approximately written as R 0a = r a N/g and R 0b = r b N/g, respectively. By taking the average over each group, we obtain the effective basic reproduction number of the infectious disease R 0 = (r a + r b )N/2g = r N/g [33] where r is the average value of the disease transmission rate of the whole population. By varying the value of r a and r b , we are able to introduce the difference in transmission rate of the infectious disease for the individuals. For the sake of comparison, we keep the average value of the transmission rate fixed as r = (r a + r b )/2. Denoting r a /r b by x, the relative disease transmission rate for the two types of individuals, after some simple algebra we have

When x is close to zero, there exists a great difference between the individuals in group A and those in group B in acquiring the disease (i.e., we consider the case where the youths are very resistant to the infection, while the adults are very vulnerable to the disease). As x goes to unity, the variation of the disease transmission rate among the two groups vanishes. Let us show in Fig. 2 the influence of the cross contact coefficient on the epidemic spreading in the population without a vaccinated compartment. In the case of the limit x → 1, we have r a ≈ r b , which means that the possibilities of acquiring the disease through susceptible-infected contact for the individuals from the two groups are almost the same. As a consequence, the final epidemic size f R , i.e., the average fraction of recovered individuals in the whole population, does not change much as the parameter varies. Note that with the current parametrization settings the final epidemic size without vaccination is about 89.3% for x = 1.0 [13] . As x diminishes, f R decreases considerably. This point can be understood by considering the case of → 0. In such a case, as demonstrated in Appendix V, due to the concavity of f R as a function of R 0 , the decrease of epidemic size f Ra in group A cannot be offset by the increase of f Rb in group B and consequently the final epidemic size of the whole system will decrease continuously as x decreases. In particular, when the value of x is less than 0.25, the value of R 0a will be smaller than unity, which means that the epidemic cannot spread throughout group A. Hence f R of the whole population is mainly contributed by f Rb and converges approximately to a value ≈ 0.5 for x < 0.25. With the increment of , the more frequent contact between the two groups will infect more individuals in group A, while the somewhat less frequent contact among those individuals from group B has just a slight impact on the final f Rb (see Appendix V). The introduction of heterogeneity of the infection rate can greatly suppress the prevalence of the infectious disease.

We now incorporate the vaccinated compartment into the epidemic spreading in the metapopulation model. We denote by f V a the proportion of the population initially vaccinated in group A. In our work we assume the same fraction of initially vaccinated individuals for the two groups, that is,

For given values of f V , x, and , we obtain the final epidemic size by implementing stochastic simulations as described in Sec. II. The simulation results are summarized in Fig. 3 , which are in good agreement with those predicted by numerically solving Eqs. (4)-(6). The overall result is that with the involvement of the vaccinated compartment, the final epidemic size will gradually decrease with the increase of f V , which is expected since vaccination can provide perfect immunity to the infectious disease and a sufficiently large fraction of vaccinated individuals can completely prohibit the propagation of the infectious disease. Though the difference between f R for x = 1.0 and that for x < 1.0 is vanishing in the limit of large f V , there exists a qualitative difference for the variation. When the individuals from the two groups interact quite frequently = 0.5, the smaller the relative disease transmission rate x is, the smaller the final epidemic size f R is. Such a dynamic scenario, however, changes when the interaction frequency among the in-dividuals from distinct groups is decreased. Specifically, a crossover behavior of f R as a function of f V emerges as the parameter drops close to zero. We notice that there arises a critical value of f V , say, f V c (whose value is about 0.45), below which the presence of heterogeneity in infection rate for the individuals from different groups can hinder the epidemic spreading, while above which the opposite effect takes place (see Appendix VI for more details). It is worth pointing out that for sufficiently small , the individuals in the two groups almost interact with others within the same group, which leads to the clustering of susceptible individuals with a high infection rate of the disease (in group B). Consequently, the disease prevalence is enlarged as compared to the case of a homogeneous interaction pattern of the two groups [e.g., the curve for x = 0.02 in the case of = 0.1 is always above that in the case of = 0.5 (not shown here)].

Now we study our model in a spatially structured population, where the individuals are located on a square lattice. For the sake of comparison, we calibrated the epidemic parameters to ensure that the infection risk in an unvaccinated population (without variation of infection) is equal across all population structures, that is, f R for x = 1 in the case of spatially structured population should be the same as f R for x = 1 in the case of a metapopulation. The simulation results are displayed in Fig. 4 , from which we note that the final epidemic size f R decreases much more rapidly as compared to that in the metapopulation case when the vaccination level increases. When the two types of individuals are randomly prearranged, f R decreases monotonically as the variation of infection increases for each vaccination level. Noticeably, we find that the crossover behavior of f R as a function of f V still exists when the interaction frequency between the two types of individuals reduces to a very low level. From Fig. 4(b) we can see clearly that there is a crossing point near f V c = 0.1. For f V < f V c , the heterogeneity in infection can efficiently hinder the disease spreading, while it promotes the propagation for f V > f V c , similar to the results in Fig. 3(b) obtained for the metapopulation model.

In what follows we investigate how the vaccination dynamics (i.e., we allow the individuals to change their vaccination behavior based on previous experience [13] ) affects the epidemic spreading in structured populations. In the initial state, we randomly choose half of the population to get vaccinated. At the end of each epidemic spreading season, we give the individuals a chance to update their strategies for vaccination before the new one starts. We implement a pairwise comparison process for the strategy updating. Specifically, whenever an individual i updates one's vaccination strategy, one just chooses an individual j randomly from one's nearest neigh- Fig. 1(b) and (b) regularly arranged as in Fig. 1(c) . The parameters are the total population size N =100×100, average value of the disease transmission rate r = 0.46 day −1 person −1 , recovery rate g = 1 3 day −1 , and number of initial infection seeds Ia=I b =10. Simulation results are averaged over 100 independent runs. bors to compare their cost (or payoff) and then adopts the vaccination choice of j with a probability dependent on the payoff difference [34] [35] [36] 

where P i and P j correspond to the payoffs of the two involved individuals and β denotes the strength of the selection. Unless otherwise specified, we select β = 1.0, implying that better-performing individuals are readily imitated, but it is not r a n d o m l y a r r a n g e d r e g u l a r l y a r r a n g e d x =0.5 r a n d o m l y a r r a n g e d r e g u l a r l y a r r a n g e d r a n d o m l y a r r a n g e d r e g u l a r l y a r r a n g e d x =0.5 We plot in Fig. 5 the epidemic size f R and the vaccination level f V in the steady state as a function of the relative cost for vaccination c for two differently arranged populations on square lattice. From Figs. 5(a)-5(c) we observe that as the value of x goes down, i.e., the heterogeneity in infection rate for the two types of individuals becomes more notable, the final epidemic level in the randomly arranged population (the open symbols) changes much more evidently than that in the case of regularly arranged population (the closed symbols). In particular, for x = 0.5, the final f R in the randomly arranged population is always greater than that in the regularly arranged population as c increases, albeit the vaccination level in the former case is slightly larger than that in the latter case for c < ∼ 0.25 [see Fig. 5(d) ]. For x = 0.3, in the randomly arranged population, though the growth trend of f R is more apparently for small c, it attains at a smaller level for large enough values of c (when the vaccination level evolves to zero), which is comparable to the case of a regularly arranged population. As x decreases even to 0.1, f R in a randomly arranged population can just grow to a much lower level as compared to that in the case of a regularly arranged population, despite the fact that the vaccination level is zero for most c values [see Fig. 5(f) ]. The reason is that the A-type individuals are difficult to infect even though they did not receive a vaccine when x is too small and as such they play the role of a natural obstructer to prevent large-scale spreading of the disease in the population. In addition, those unvaccinated A-type individuals will attract other individuals to not get vaccinated, giving rising to very low level of vaccination in the steady state [Figs. 5(e) and 5(f)]. For a regularly arranged population, however, since the B-type individuals are clustered together, they are very prone to the attack of disease, and consequently the final epidemic can reach a rather large level.

In summary, we have incorporated the heterogeneity in infection rate of individuals and also the vaccination dynamics into the traditional susceptible-infected-recovered compartmental epidemic model to study their potential effects on the disease prevalence and vaccination coverage. For this pur-pose, we have considered a more practical framework where the whole population is classed into two types of groups whose members are endowed with different capabilities in catching a disease. To keep things simple, the individuals within the same group are assumed to be identical in their infection rate. The proposed model has been investigated in a simple metapopulation and spatially structured populations, with and without involvement of vaccination, by using numerical simulations as well as analytical treatments.

We have shown that whether the introduction of heterogeneity in the infection rate of the individuals exerts positive or negative effects (i.e., hampers or expedites) on the epidemic spreading depends closely on both the extent of the heterogeneity of the disease transmission rate and the interaction frequency among the individuals from different groups. To be more specific, the heterogeneity in infection rate can always give rise to a decrease of the final epidemic size provided the individuals from different groups interact with equal likelihood. Nonetheless, as the individuals become more inclined to interact mainly with others from the same group, the heterogeneity in infection rate can hinder the epidemic spreading only in the situation that the fraction of individuals vaccinated is low enough. Very surprising, this just facilitates the epidemic spreading in a regime with the presence of a large fraction of vaccinated individuals (but not large enough to eradicate the disease completely).

Our work is expected to provide some valuable instructions for the prediction and intervention of epidemic spreading in the real world. The results summarized in Figs. 2-5 suggest that when evaluating the seriousness of an epidemic, we should take into account both the factors of the diversity of the infection rate of the individuals and the interaction patterns among them simultaneously; otherwise we may overestimate or underestimate the spreading trend. Alternatively, without such considerations, we may overshoot or undershoot the desired amount of action when developing, regulating, and making vaccine policy. In addition, when individuals are allowed to change their vaccination decisions according to their experience and observations, we find that as the heterogeneity in infection rate for the two types of individuals becomes more noticeable, the final epidemic level in randomly arranged population changes much more evidently than that in the case of a regularly arranged population, hence giving us a vital clue as to how to make efficient vaccine campaign, namely, we should distribute the vaccine in the population as widely as possible so that the spreading path of the disease can be efficiently suppressed.

To summarize, our proposed model captures essential elements in real-world epidemic spreading, which has not been fully discussed previously. Therefore, we believe our results will give some insights to the policy makers. There are still many issues, such as diversity of recovery rate, heterogeneous cost for infection and vaccination, and more complex contactnetwork structures, which are totally overlooked in the present work and deserve to be explored in the future. In addition, the spread of awareness of the epidemic and/or the vaccination sentiment would also impact greatly the vaccination behavior of the individuals and hence the epidemic outbreaks [37] [38] [39] .

We hope our work could stimulate further work in this line of research.

This work was supported by the National Natural Science Foundation of China (Grants No. 11005051, No. 11005052, and No. 11135001 ).

Here we present theoretical analysis for the simple metapopulation model. For convenience, let us denote r N/g by C, which is kept as a constant. The combination of Eqs. (1)-(3) with Eq. (7) yields

After eliminating Ia I b and I b Ia from these equations we readily obtain

Now we integrate these two equations with respect to time from 0 to ∞. By using the initial condition S a (0) = S b (0) ≈ 1 and R a (0) = R b (0) = 0 and the final state I a (∞) = I b (∞) = 0, we get the following two transcendental equations for the final epidemic size R a (∞) and R b (∞) for each group:

] .(16) What we want to figure out is the relationship between the final epidemic size f R and the cross coefficient , so we take a derivative of Eqs. (15) and (16) with respect to and get

After doing some algebra we obtain We can rewrite this equation as

where

In the case of = 0, it is easy to verify numerically that K > 1 for all our x values of interest (say, x > 0.01). More intuitively, for SIR model in a well-mixed population, the final epidemic size is determined by the self-consistent equation R(∞) = 1 − exp −R0R(∞) . Figure 6 features the solutions, from which we note that f R is a concave function of R 0 . If we decrease the value of x such that in the limit of = 0 the variables R 0a and R 0b always satisfy the relationships R 0a < R 0b and (R 0a + R 0b ) /2 = C = R 0 , then due to the concave curvature, the variation of the final epidemic size in group A will be more remarkable than that in group B, as illustrated in Fig. 6 . For each fixed value of x, as increases, the increasingly frequent contact among the individuals from different groups will have a greater affect on R a (∞) than on R b (∞) (as long as x is not too small), giving rise to the increase of f R in the whole population. Since R a (∞) increases and R b (∞) decreases with the increment of , the value of K will decrease monotonically according to Eq. (21) , which is reflected correctly in Fig. 2 .

Here we demonstrate the existence of the crossover behavior for the curves of f R as a function of f V = y for different values of x. In a well-mixed population, we know the final fraction of recovered population for the SIR model satisfying the equation R(∞) = 1 − exp −R0R(∞) . When (5) a Results of (f V c , f R ) obtained from stochastic simulations. b Results of (f V c , f R ) predicted by analytical treatments. a proportion y of preemptive vaccination in introduced before the epidemic starts, we can readily obtain R(∞) = (1−y) 1 − exp −R0R(∞) . For our proposed model, we consider two limited cases. The first case is x = 1, i.e., the individuals in the two groups are identical, and in such a case we have

where R 0b = r b N/g = r N/g = C.

The other limited case is x → 0, which means that the disease transmission rate for the individuals in group A is nearly zero. By approximating R a (∞) | x→0 = 0 and combining Eqs. (15) with (16) we have

where R 0b = r b N/g = 2 r N/g = 2C. We assume that the curves of f R for the two cases have a crossing point so that

Denoting R b (∞) | x=1 by z, combining Eqs. (22) , (24) , and (26) , and recalling that C = 2.5, we obtain exp −10(1− )z = 2 exp −2.5z −1.

To validate the assumption, Eq. (27) must have an exact solution, which means that

Solving the inequality yields < 0.5. That is to say, the crossover behavior will always exist as long as is strictly smaller than one-half. From Eqs. (27) and (22) 

dy dz = exp −2.5z +2.5z exp −2.5z −1 (exp −2.5z −1) 2 < 0.

By dividing Eq. (30) by (29) we get d /dy > 0, which indicates that the crossing point will move to the right (i.e., the curves intersect at larger values of y = f V ) with an increase of the cross contact coefficient . In Table I we summarize the crossing point values (f V c , f R ) of the curves for x = 1.0 and 0.02 yielded by the stochastic simulations as well as those predicted by Eqs. (27) , (22) , and (24) . We notice that the results obtained from different methods match quite well with each other. The invisible differences may be due to the finitesystem-size effect. Specifically, with increasing the curves for x = 1.0 and 0.02 intersect at points with larger f V c .

Proc. Natl. Acad. Sci. USA

Infectious Diseases of Humans: Dynamics and Control

Modeling Infectious Diseases in Humans and Animals

Note that only in the limit case of → 0 can R0 be approximately written as (ra + r b )N/2g. In any cases 0, we are unable to write out the explicit form of R0

Proc. Natl. Acad. Sci. USA

Proc. Natl. Acad. Sci. USA