Assessing Intervention Strategies for Non Homogeneous Populations Using a Closed Form Formula for R0
AAssessing Intervention Strategies for NonHomogeneous Populations Using a Closed FormFormula for R • Wolfgang H¨ormann • Refik G¨ull¨u
Bogazici University Industrial Engineering Department, Bebek 34342 Istanbul [email protected] • [email protected] • refi[email protected] Abstract
A general stochastic model for susceptible → infective → recovered (SIR) epidemics innon homogeneous populations is considered. The heterogeneity is a very important aspecthere since it allows more realistic but also more complex models. The basic reproductionnumber R , an indication of the probability of an outbreak for homogeneous populationsdoes not indicate the probability of an outbreak for non homogeneous models anymore,because it changes with the initially infected case. Therefore, we use “individual R ”that is the expected number of secondary cases for a given initially infected individual.Thus, the effectiveness of intervention strategies can be assessed by their capability toreduce individual R values. Also an intelligent vaccination plan for fully heterogeneouspopulations is proposed. It is based on the recursive calculation of individual R values. Keywords: basic reproduction number, stochastic epidemics, susceptible infected re-covered, non homogeneous populations, intervention methods, intelligent vaccination.
Epidemiological information is mostly used to plan and evaluate strategies that preventdisease spread by identifying risk factors. Therefore, various disease spread models weredeveloped in the literature. This paper considers a stochastic model for susceptible → infective → recovered (SIR) epidemics among non homogeneous populations. The basicreproduction number R as the expected number of secondary cases produced by a singleinfected case in a totally susceptible population is important for the determination of theoutbreak probability under homogeneous mixing assumption (Allen and Burgin, 2000;Craft et al. , 2013; Hernandez-Suarez, 2002; Kumar and Goel, 2020). In order to assessthe dynamics of disease behavior, Ross (2011) analyses the effect of population size and1 a r X i v : . [ q - b i o . P E ] A ug isease time distribution on R . Moreover, heterogeneity in population has importanteffects on disease spread behavior as discussed eg. in (Inaba, 2012; Meyers, 2007) withits realistic assumptions.This paper is concerned with the notion of computable R for heterogeneous models.Our aim is to calculate the expected number of secondary cases produced by a uniquegiven infected case and use it to develop effective intervention strategies and assess theintervention strategies without simulation. Watson (1980) proposes a model with twomixing levels and defines three types of outbreaks: localized, restricted and generalizedoutbreaks. There are also others considering R for non homogeneous populations toestimate outbreak probabilities. Since the definition of R for heterogeneous populationsis complex, Diekmann et al. (1990) and Trapman (2007) estimates bounds for the ex-pected infectivity based on disease length and heterogeneity state regarding deterministicdisease spread models among totally heterogeneous populations. Because the analysisof stochastic disease spread models for heterogeneous populations is difficult, either sim-ulation is used for the analysis (Ajelli et al. , 2010; Longini et al. , 2005) or the modelsare simplified by decreasing the number of mixing levels. There are a number of recentpapers on epidemics in heterogeneous populations in order to estimate R . They imple-ment agent based simulation claiming that R cannot be directly calculated (Lipsitch etal. , 2003; Longini et al. , 2004) while Ball (1986) calculate the average R exactly amonga stratified population with two levels by applying branching process methods. Keeganand Dushoff (2016) estimate R for finite populations with different heterogeneity typesto analyse the effects of heterogeneity on basic reproduction numbers. However, a single R value is calculated in all of these models to analyse the possibility of an outbreak.Artalejo and Lopez-Herrero (2013) measures R e, that is the exact number of secondarycases generated by the tagged infected individual if the epidemic starts or is already inprogress. They emphasize the probability distribution of the number of secondary casesrather than its mean. Following Artalejo and Lopez-Herrero (2013), Economou et al. (2015) determine the distribution of the number of secondary cases for SIS models withexponential disease time. L´opez-Garc´ıa (2016) introduces R exactx as the random variabledenoting the number of individuals directly infected by a given infected individual dur-ing his infectious period given the current state of the process x and its probability mass2unction is investigated.Moreover, R is generally used in the literature for analyzing the possibility of anoutbreak even if it is also possible to use it for intervention strategy analysis. In theliterature, to develop and analyze epidemic control strategies, different mathematicalapproaches are implemented like introducing contact network epidemiology (Dimitrovand Meyers, 2010), implementing optimal control tools (Sharomi and Malik, 2017) andsimulating scenarios (Carvalho et al. , 2019; Wu et al. , 2006). However, there are somerecent papers considering the use of R to analyze and develop epidemic control meth-ods. Ball (1986) use the average R for optimal vaccination policies in a populationpartitioned into households. Artalejo and Lopez-Herrero (2013) also suggest to use R e and R p to design control strategies for prevention of an outbreak. They consider Markovchains while modelling disease spread, so they assume exponential infectious period. R e indicates the exact number of secondary cases produced by a single infective while R p de-notes the exact number of secondary cases produced by all currently infected individualsuntil first recovery. Markov chains with exponential infectious period and homogeneousmixing assumption allow them to obtain exact measures and valuable insight rather thanan expectation. Since we consider a totally non homogeneous population and a discreteinfectious period distribution, we can only measure the expected number of secondarycases produced by a unique given infected case that changes with the selected initiallyinfected individual and we call it individuals R .In this paper, we make three main contributions. Firstly, we introduce individual R as the expected number of secondary cases produced by a unique given initially infectedindividual. Individual R is an expectation rather than the exact number of secondarycases so we can propose a general formula for individual R in this paper that is applicableto all types of heterogeneous populations with any size. Our second contribution is topresent that it is possible to assess intervention strategies by using the exact formulafor individual R without reverting to simulation. It is possible to assess the impactof intervention strategies by their capability to reduce individual R values. Also, amaximum individual R value smaller than 1 guarantees that an outbreak is impossible.Lastly, intelligent intervention strategies can be identified based on individual R values.We propose a vaccination strategy such that the individual with greatest individual R are3accinated first. In order to choose the individual who is vaccinated next, we recalculateindividual R values for the unvaccinated individuals and choose the individual with thegreatest individual R again. Thus, our vaccination strategy is to vaccinate individualsone by one by choosing the susceptible having the largest individual R .The paper is structured as follows. In Section 2, the notion of R for non homoge-neous models is explained and the general formula to calculate individual R values fornon homogeneous populations is presented. In Section 3, basic non homogeneous modelsfor influenza spread are presented. Then, in Section 4, we evaluate some interventionmethods applied for these models using their individual R values. Moreover, we proposeintelligent vaccination plans based on individual R . We also consider a model with over-lapping mixing groups and check how individual R values changed due to interventionstrategies in Section 5. Finally, concluding remarks are given in Section 6. R for Non Homogeneous Models In this paper, a stochastic SIR model is considered in a non homogeneous population.Stochastic SIR models in large non homogeneous populations grew popular among prac-titioners in recent years, see eg. (Ajelli et al. , 2010; Longini et al. , 2005). The infectionprobability between an infected and a susceptible individual is modelled with a compar-atively small number of parameters assuming mixing in overlapping mixing groups. Thedetailed structure of a population is generated such that the mixing groups match in sizeand age those of real world census data. As mixing groups typically households, neigh-bourhoods, communities, schools and work places are considered. In several papers it isassumed without any discussion that the only way to assess the behaviour of such modelsis simulation. This fact attracted our attention and we aim to develop here an approachto assess the behaviour of such models for large populations using a properly definedbasic reproduction number R that can be calculated easily also for large populations.As individuals in a non homogeneous population are not identical, R for non homo-geneous populations depends on the initially infected individual that is chosen. Thus,different R values occur for different initially infected individuals. In agent based simu-lation literature, the value of R for the entire non homogeneous population is generally4stimated by assuming “the random case”. Thus, the initially infected individual is se-lected among the population with equal probability for every individual. Then, they usean average to calculate R (Longini et al. , 2005). It is also suggested to estimate age de-pendent R values and calculate overall R as a weighted average of age dependent attackrate patterns (Germann et al. , 2006). The studies which use branching process methodsto calculate R also considers R for non homogeneous populations as a mean of differentsecondary cases for different initially infected individuals (Ball and Lyne, 2002). In thesestudies, populations with different mixing levels and moderate size are considered basedon census data. However, average R is estimated via simulation without exact solutionand it cannot be used to assess the possibility of an outbreak anymore.That this is nota sensible approach can be demonstrated with a very simple example:A population with N = 200 is composed of two sub-populations of equal size A and B . An infective from A infects a susceptible from A with probability 0.003 and a sus-ceptible from B with probability 0.0005 during his total infectious period. Furthermore,an infective from B infects a susceptible from B with probability 0.015 and a suscep-tible from A with probability 0.0005 during his total infectious period. We can easilycalculate: The expected number of secondary cases for a single starting infective of A is 99(0.003)+100(0.0005)=0.347 and for a starting infective of B it is 1.535. Taking theaverage over all individual we get 0.941. A value of R below one should indicate thatan outbreak is impossible, but in our little example it is clear that an outbreak in group B is likely if the first infective is of group B . And in such a case also several individualsof group A are likely to be infected. R for non homogeneous models is studied especially by using Markov models sincethey allow to calculate R exactly in the literature. However, there are some problemswith Markov modelling of disease spread. Markov chain processes requires exponentialdisease time which is clearly unrealistic. Meanwhile, the complexity of Markov modelsfor non homogeneous populations increases exponentially due to the size of state space,so the exact distribution of R ( i ) can only be calculated for very small populations. It isclear that even for moderate sized populations the state space for such a model is huge.This makes numerical calculations so difficult that L´opez-Garc´ıa (2016), who considers asimilar but continuous time model with exponential disease times and develops numerical5ethods to calculate the distribution of important stochastic descriptors, stresses evenin the title of the paper that this is only possible for small networks.Following the simulation literature, we consider a simple discrete time stochasticmodel. Important is that we allow a very general mixing structure assuming that ina finite population of size N we know all probabilities p ij that within one time-step (inpractice typically one day) an infected individual “ i ” transmits the disease to a susceptibleindividual “ j ”. It turns out that it is also sometimes necessary to allow the possibilitythat the infection probabilities change with time. In such cases we will write p ijt .The estimation p ij for each pair of i and j is possible for small population sizes likehospitals etc. Laskowski et al. (2011) implement an agent based modelling for the spreadof influenza like disease in an emergency department. They model patients as occupyinga circular space with a radius of 60cm and define different contact types like close andcasual contacts based on the distance between agents. Moreover, they consider a basicpatient flow model throughout which agents come into contact with each other andthe probability of infection is found based on the agent distance during contact and theduration of the contact. However, the estimation of p ij is very difficult for large populationsizes. The overlapping mixing groups approach is mainly suggested for estimating p ij inlarge populations. Individual based models for disease spread have been implementedduring the last 50 years, but it has been popular recently due to the lack of both dataand advanced computational availability in the past Yang et al. (2008). Carley et al. (2006) propose a scalable city wide multiagent network numerical model where agentsare embedded in social, health, and professional networks. The model allows to defineheterogeneous population mixing by agent and social networks characteristics based onreal data from census, school districts, general social surveys, etc. Bian (2004) presentsa conceptual framework for an individual based spatially explicit epidemiological modelbased on the following assumptions: (1) individuals are different so age groups are needed;(2) an individual has contacts with a finite number of individuals in different clusters likehome and workplace; (3) individuals travel between clusters; and (4) the individualshave different contact rates such as fewer contacts for retired individuals than employedindividuals. Thus, two types of contacts are defined: those within a group and thosebetween groups. Moreover, the shift from population based models to individual based6odels is explained by the rapid improvements in computing power and availability ofspatial data. Longini et al. (2004) compare the efficiency of the use of anti viral drugs andvaccination for a population with 2,000 persons who are stochastically generated by theage distribution and approximate household size published by the US Census Bureau.Yang et al. (2008) study an individual space time activity model for the target city,Eemnes in the Netherlands based on an activity survey data, a synthesized householddata, land use data, and PC6 statistical data.The agent based models collect the infection data at individual level and becomemore realistic. However, its increased complexity also brings too much a burden formodel structure and requires simulation. Moreover, p ij s are required to be computedfor individual based models by considering the infection probabilities p f , p s , p n , and p c in mixing groups of “households”, “school and play groups”, “neighbourhoods” and“communities”, respectively that are also changing with age groups. Then, if infectionevents between different mixing groups are assumed to be independent, the probabilityof infection between two individuals i and j during a day is calculated as p ij = 1 − (1 − p c ) I C ( j ) (1 − p n ) I N ( j ) (1 − p w ) I W ( j ) (1 − p f ) I F ( j ) (1)where the indicator function of a subset is defined as I M ( j ) = (cid:40) j ∈ M ixing Group M of i otherwise. That implies that individuals i and j can mix in different mixing groups in set M (com-munity, neighborhood, school, work, family etc.). So the indicator function is one forseveral j.The state of our model is described by the state vector holding the state S , I or R for all N individuals. In one time step a susceptible individual j is infected by a singleinfected individual i with probability p ij . If there is more than one infected individualthe assumption that these infections are independent of each other leads to the totalinfection probability for individual j : p j = 1 − (cid:89) i ∈ I (1 − p ij ) , I denotes the set of infected entities. The new infections are thus a sequence of | S | independent Bernoulli trials with probabilities { p j | j ∈ S } , where S denotes the set of allsusceptible individuals. To pass from state I to R we use the model assumption that thedisease times of all individuals are independent and follow a discrete distribution withprobability mass function f D ( d ) with d = 1 , , . . . .For non-homogeneous mixing it is obvious that we need, like suggested in L´opez-Garc´ıa (2016), a definition of R that considers which individual is the single startinginfected. As we consider here large populations, we use the simple classical definition of R and define: R ( i ) = E[secondary cases for starting with a unique infected individual i ]and call it individual R .One important advantage of individual R is that it can be calculated easily also forlarge populations. To develop the formula we first have to calculate the probability (cid:101) p ij that susceptible j is infected by infectious i during the total disease time of i . This iseasily done using conditioning on the disease time D : (cid:101) p ij = ∞ (cid:88) d =0 [ f D ( d )(1 − (1 − p ij ) d )] . (2)It is also sometimes possible that the infection probabilities change with time written as p ijt . Then, the probability (cid:101) p ij can be calculated as (cid:101) p ij = ∞ (cid:88) d =1 [ f D ( d )(1 − Π dt =1 (1 − p ijt ))] . (3)Note that also for a disease time distribution with unbounded domain it is not difficultto calculate a close approximation of (cid:101) p ij as the error commited by a cut off of the sum after d = d m is obviously always smaller than 1 − F D ( d m ) and can thus be easily controlled.Individual R ( i ) is then simply the ”column sum of the matrix (cid:101) p ij ” or more precisely thesum of all (cid:101) p ij ’s for i fixed and j = 1 , , . . . , i − , i + 1 , i + 2 , . . . , N : R ( i ) = (cid:88) j : j (cid:54) = i (cid:101) p ij . (4)The complexity of calculating R ( i ) in (4) for i = 1 , , . . . , n is in total O ( d m N ),where d m denotes the size of the domain of the disease time D for bounded disease timeor the cut off value of the infinite sum for the case that D has an unbounded domain.8 .1 Use of Individual R on Intervention Analysis A main aim of building agent based simulation models for influenza spread is the assess-ment of interventions. How is the spread of the disease changed for instance, when •
15 percent of all individuals are vaccinated; • anti-viral drugs are given to all members of a household when one member turnsout to be infected; • when 50 percent of all infected would stay at home after the first day of the disease.How can the calculation of all R ( i ) values help to asses the behavior of the diseasespread? As we have demonstrated with the help of a simple example above, the average ofall R ( i ) values does not allow a direct assessment. But it is easy to see that if max i R ( i )is smaller than one, an outbreak is impossible. If that value is above one the behavior isnot certain but an outbreak is possible.Like for many other interventions also for the second and third intervention exampleabove it is obviously necessary to assume that the p ij values change with time and aredenoted by p ijt on day t . The (cid:101) p ij ’s are obviously calculated using Equation 3. To obtainthe R ( i ) we need again the column sums given in (4).To quantify the influence of such interventions it is first necessary to decide how theparameters of the model are changed by the intervention. Here it may be necessary tomake assumptions (or guesses) how the infection probabilities are changed; if we considerthe case that when 50 percent of all infected would stay at home after their first day ofinfection is an example where it is clear that people staying at home have infectionprobabilities of zero with all individuals not belonging to their household. It is possible to calculate individual R values exactly for all non homogeneous modelsusing Equation 4. In this part, we calculate (cid:101) p ij and individual R values for some nonhomogeneous population structures in the literature that are well applicable for airborne9iseases like influenza. We need a discrete disease time for influenza and assume like ? that the probability mass function of 3, 4, 5 and 6 days with probabilities 0.3, 0.4, 0.2,and 0.1 respectively . Moreover, we consider two different non homogeneous populationmodels. Then, in Section 4, we evaluate some intervention methods applied for them. We consider a network of cities around the world connected by transportation. Thismodel is commonly referred as meta population model in the literature suggested byLevins (1968). It includes several sub populations in which perfect mixing is assumed.Individuals travel between the cities leading to disease spread according to probabilisticrules based on the population size and the travel frequency between the cities. Populationsize and travel data can be obtained from different available sources (e.g. PopulationDivision, U.S. Census Bureau 2004).Consider now three cities, numbered 1, 2 and 3. Assuming symmetry in travel, weconsider a function p ij given in Table 1 and compute individual R values by applyingEquation 4. In big cities, it is standard to assume that R is the same in a homoge-nous population. Thus, the expected number of individuals infected by a single infectedindividual is not increasing with the size of the population and the probability to meetand potentially to infect another individual is reduced as the population size increases(Lund et al., 2013). Therefore, we assume greater infection probabilities within a citywith smaller population sizes. To obtain the infection probabilities between cities is morecomplicated and challenging and it is beyond the scope of this work (Lund et al. , 2013).Here, we take a simplistic view and assume that the travel frequency is the greatestbetween city 1 and city 2 and the smallest between city 1 and city 3 by considering thedistances between cities and population sizes. Furthermore, the infection probabilitiesbetween cities are considered to be around 2 percent of the infection probabilities withinthe same city. However, it is also possible to obtain travel frequency data for better es-timation. Moreover, the reported R values for the basic reproduction number in a fullysusceptible population is in the range of 1.6 to 2.4 for influenza (Germann et al. , 2006).Thus, while setting the infection probabilities, we target to obtain average R et al. (2004). We estimate the infection probabilities by dividing10 op.Size City 1 City 2 City 3 SuperSpreaders R ( i )City 1 746 5.20e-4 1.30e-5 8.66e-6 1.00e-3 1.673City 2 500 1.30e-5 7.80e-4 1.04e-5 1.50e-3 1.714City 3 746 8.66e-6 1.04e-5 5.2e-4 1.00e-3 1.668Super-spreaders 8 1.00e-3 1.50e-3 1.00e-3 0 9.174 Table 1: Population Sizes and Infection Probabilities for The Population with MultipleCities.target R = 1 . R values in this case is equal to the number of cities plusone for the super spreaders. The corresponding individual R values are given in the lastcolumn of Table 1 computed by using Equation 4 consistent with the simulation results. We also consider a population partitioned into several households similar to Ball and Lyne(2002) since the household based public health interventions are important to preventthe spread of infectious diseases. Moreover, the two levels of mixing is also important forthe behaviour of the epidemic.Lets consider that an infected individual infects a household member with probability p h and other individuals with probability p c . p h is selected considerably higher than p c since individuals in the same household have closer contacts. If we denote the familymembers of individual i as set N ( i ), the infection probabilities for individual i are p ij = (cid:40) p h , if j ∈ N(i) p c , otherwise . (5)We assume that p h and p c are 0.0001 and 0.06 respectively for our intervention analysisand we consider 498 households each consisting of four individuals. Further, there mightbe some individuals in the population who meet with other people more frequently thanother individuals eg. due to their work. We call such people super spreaders. In theliterature, super spreaders are defined as the individuals infecting more contacts than11thers. We assume that each infected super spreader infects with probability p s = 0 . Intervention methods aim to change the characteristics of the spread of a disease bychanging the infection probabilities p ij . We suggest to assess the impact of interventionstrategies by calculating and comparing the individual R values of the different scenarios.We consider the models described in section 3 where the individuals within the samegroup are assumed to behave homogeneously. Colizza and Vespignani (2008) define theusual R as a function of disease parameters for each group while a subpopulationsreproductive number R ∗ as a function depending on the diffusion rate of individualsamong subpopulations. Thus, a group specific basic reproductive number is consideredfor a deterministic metapopulation system and the epidemic behaviour on both the globalscale and the local scale is determined by R ∗ and R , respectively. Barth´elemy et al. (2010) consider a stochastic metapopulation model by taking into account both temporaland topological fluctuations. Moreover, individual R computed in this section is alsoa group specific basic reproduction number by considering both infection among thepopulation members of each group and between the members of different groups instead oftwo different basic reproduction numbers as in the study of Colizza and Vespignani (2008).However, individual R values can be generalized for every non homogeneous populationmodel like individual based models. Furthermore, we illustrate some numerical resultsto demonstrate the use of individual R for both developing and assessing interventionstrategies including vaccination, social distancing and use of antiviral drugs. For the evaluation of vaccine efficacy, it is assumed that vaccination takes place be-fore the infection starts to spread and that all vaccinated individuals develop immunity.Therefore, vaccinated individuals are not considered as susceptible anymore. For thevaccination as an intervention strategy, it is possible to assume random vaccination inwhich the individuals who are vaccinated are selected randomly with equal probabilities12ithin the population. However, it is better to use the vaccine efficiently to attain herdimmunity by vaccinating a smaller number of individuals.Ball and Lyne (2002) develop optimal vaccination policies for a population with twolevels of mixing and consider optimality in terms of the cost of vaccination programincluding vaccine, administration, and travel. Here, we propose an intelligent vaccinationstrategy when assuming that the cost of vaccine is considerably larger than the cost ofvaccination. In other words, the aim is to obtain for a fixed number of vaccines thegreatest reduction for the maximum individual R value . In this vaccination strategy,individuals with large individual R are vaccinated first because we both eliminate thegreatest individual R and obtain the greatest total reduction in the other individual R values if p ij s are symmetric. Therefore, as a next step all individual R values must berecalculated and their values are arranged in non increasing order. Then, the individualwho is vaccinated is selected from the top of the list and the individual R values forunvaccinated individuals are recalculated. Thus, our intelligent vaccination policy is tovaccinate individuals one by one choosing the susceptible having the largest individual R .By taking the population matrix, popm and the target number of vaccinated individuals, v critial as input parameters, Algorithm 1 presents the intelligent vaccination strategy. Algorithm 1
Intelligent Vaccination Strategy Set v=0 Compute ˜ p ij using Equation 2 for i = 1 , , . . . , N − v do Compute R ( i ) using Equation 4 end for Order R ( i ) from largest to smallest Remove individual i with the largest R ( i ) from the population matrix Set v = v + 1. If v < v critical go to step 3. Otherwise, stop the algorithm.In the theory of branching process where m is the expected number of children of eachindividual, m < m values are different for different individuals(Antreya, 2006). If the maximum m is smaller than one, then the process will be alsoextinct with probability one. Since we consider non homogeneous populations yieldingdifferent individual R values, we guarantee that there will be no outbreak by reducing13ll individual R values below one. Algorithm 1 for the intelligent vaccination policyis a greedy heuristic for heterogeneous populations and it approximates to the optimalvaccination policy as the heterogeneity level decreases.If we consider the population with three cities described in Section 3.1, the intelligentvaccination strategy requires to vaccinate individuals from different cities. The simula-tion results indicate that a significant proportion of the population has been infected withprobability 0.662 without vaccination. To guarantee that the infection is going to disap-pear before involving a significant number of the population by implementing Algorithm1, we observe that individual R values in all cities reduce below 1 if 292 individuals fromcity 1, 200 individuals from city 2, 290 individuals from city 3 and all 8 super spread-ers are vaccinated. Therefore, the minimum number of required vaccinated individualsreducing all individual R values to under 1 is found to be 790 where there will be nooutbreak controlled by computing the final outbreak size through simulation. As ’in thecity infection probabilities’, p ii are considerably greater than ’between the cities infectionprobabilities’, p ij where i (cid:54) = j for a model with multiple cities, intelligent vaccinationstrategy based on sequential vaccination also gives us the optimal vaccination strategyfor reducing all individual R values to under one. Let 291 individuals be vaccinated fromcity 1 instead of 292 individuals, then more than one individual have to be vaccinatedfrom the other cities to decrease individual R value of city 1 below one since p is con-siderably greater than p and p . This also holds for city 2 and city 3. However, it doesnot always yield the optimal strategy. If we consider an individual based model whereeach individual has its unique R ( i ), we need to decide which individuals are vaccinatedin one step by considering all relationships. Even if it is not possible to vaccinate enoughpeople to reach herd immunity intelligent vaccination is still important in order to havethe greatest possible reduction of the individual R values. Table 2 shows how individ-ual R values change for an increasing number of vaccinated individuals when using theintelligent vaccination strategy.For the population partitioned into households described in Section 3.2, the individual R value for the 1992 individuals living in households is 1.509. The sequence of intelligentvaccination starts with the super spreaders. Then, one individual is vaccinated from everyfamily. To reduce the maximum individual R value below 1, vaccination of one individual14 VaccinatedIndividuals 108 VaccinatedIndividuals 208 VaccinatedIndividuals 308 VaccinatedIndividuals 790 VaccinatedIndividualsCity Vacc.Ind. R ( i ) Vacc.Ind. R ( i ) Vacc.Ind. R ( i ) Vacc.Ind. R ( i ) Vacc.Ind. R ( i )1 0 1.644 35 1.563 73 1.479 111 1.397 292 0.9992 0 1.671 32 1.563 56 1.479 81 1.396 200 0.9993 0 1.639 33 1.562 71 1.480 108 1.397 290 0.999 Table 2: Individual R Based Intelligent Vaccination with Different Number of VaccinatedIndividuals for The Population with Multiple Citiesfrom every family is not sufficient so the vaccination continues with the vaccinationof second individuals from each family. It is easy to calculate that the maximum ofindividual R values drops below one when two individuals are vaccinated in 139 familieswhile only one individual is vaccinated from the remaining 359 families. The two resultingindividual R values are 0.999 (for 1077 individuals) and 0.777 (for 278 individuals). Theminimum number of vaccinated individuals required for herd immunity is thus foundto be 645 and. Furthermore, Table 3 presents the number of susceptibles with theircorresponding individual R values if 8, 257, 506 and 755 individuals are vaccinatedrespectively. Table 3 indicates that individual R is 1.483 for unvaccinated individualsif 8 individuals are vaccinated while individual R is reduced to 1.159 and 1.381 for747 and 996 unvaccinated individuals, respectively if 257 individuals are vaccinated. If506 individuals are vaccinated, individual R is reduced to 1.057 for all unvaccinatedindividuals. Moreover, the last two columns of Table 3 show that individual R s arereduced to 0.732 and 0.995 for 498 and 747 unvaccinated individuals if 755 individualsare vaccinated, so vaccinating more than 645 individuals decreases individual R muchlower than 1.If intelligent vaccination is compared to random vaccination, it is observed that theminimum number of required individuals to be vaccinated to reach herd immunity is muchhigher for random vaccination and its performance under limited vaccination supply isalso clearly worse. 15 IndividualsVaccinated In 50% familiesone membervaccinated In all familiesone membervaccinated In 50% families twoand in 50% families onemember vaccinated R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv.0 8 0 257 0 506 0 7551.483 1992 1.159 747 1.057 1494 0.732 498- - 1.381 996 - - 0.955 747 Table 3: Individual R Based Intelligent Vaccination with Different Number of VaccinatedIndividuals for The Population Partitioned into Households
WithoutQuarantine Quarantine withCompliance Rate 50% Quarantine withCompliance Rate 80% R ( i ) Number of Ind. R ( i ) Number of Ind. R ( i ) Number of Ind.1.509 1992 1.191 1992 0.999 19928.084 8 0 8 0 8 Table 4: Quarantine After First Day of Infection for The Population Partitioned intoHouseholds
The simplest intervention strategy that can be considered as a method of social distancingis household quarantine. The effectiveness of household quarantine depends on manyadditional disease parameters like the time between the start of the infection and thestart of the symptoms and the compliance rate indicating the percentage of symptomaticinfluenza cases who remain at home. Household quarantine can be implemented onlysome time after the infection starts so we assume that it is implemented after the firstday of the disease.We consider the population partitioned into households only since it is not possible toimplement social distancing by the nature of a model with multiple cities. To demonstratethe impact of household quarantine, it is assumed that individuals stay at home after thefirst day of infection with probability 0.5 suggested in the study of Wu et al. (2006). Table4 shows the resulting changed individual R values. The important point in Table 4 is thereduction in the individual R values of the household members. Therefore, it is possibleto decrease individual R values by increasing compliance rate. It may be possible to16 R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv.1.448 1992 1.386 1992 1.323 1992 1.258 19927.942 8 7.797 8 7.649 8 7.498 8 Table 5: Use of Anti Viral Drugs with Different Reduction Factors without HouseholdQuarantine for The Population Partitioned into Householdsincrease the compliance rate if a viable diagnostic support including virological testing isavailable. Thus, we search the compliance rate to attain herd immunity for the householdmodel. We observe that household quarantine must be accepted by at least 80% of theinfected to make an outbreak impossible for the population partitioned into householdsdescribed in Section 3.2.
Antiviral drugs can be both of prophylactic and therapeutic importance. The use ofantiviral drugs that is evaluated here prevents infection given exposure. Therefore, it isassumed that antiviral drugs reduce the probability of transmission to others and theprobability of being infected given exposure. There are no direct estimates of how muchantiviral drug will reduce the probability that an infected individual will develop in-fluenza symptoms compared with an infected person who is not using antiviral drugs butthese parameters are inferred from household studies of antiviral drugs in the literature(Longini et al. , 2004). Therefore, considering that family members of the initially in-fected individual use antiviral drugs, we check how their individual R values change forassuming different reduction factors of anti viral drugs. The results in Table 5 indicatethat, as expected, the effectiveness of the use of anti viral drugs is strongly influencingto the reduction capability for infection probabilities.Furthermore, we check how effective is the combination of anti viral drugs and house-hold quarantine. The results are given in Table 6. We observe that assuming a compliancerate of 50% and reduction rate 40% it is possible to prevent an outbreak by using thecombined strategy even if this is not possible when using household quarantine and anti17 R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv. R ( i ) Num. ofIndiv.1.130 1992 1.068 1992 1.005 1992 0.940 19925.477 8 5.333 8 5.184 8 5.034 8 Table 6: Use of Anti Viral Drugs and 50% Household Quarantine with Different Reduc-tion Factors for The Population Partitioned into Householdsviral drugs alone.
The overlapping mixing group model tries to imitate the disease spread in a real worldcommunity using census data. It requires only a moderate number of parameters. Theaverage R for these models is calculated using simulation in the literature (see Longini et al. (2004)). The model uses several mixing groups like “households”, “school and playgroups”, “neighborhoods” and “communities” with their respective infection probabilities p f , p s , p n , and p c changing with age groups to model all infection probabilities p ij thatcan be calculated by using Equation 1.Moreover, following Longini et al. (2004) we also consider asymptomatic cases for theoverlapping mixing group case as a feature of influenza in real world. Asymptomatic casesare the infected individuals who do not have symptoms. Their infection probabilities arealso considered to be smaller than the ones for symptomatic cases. The implementationof intervention strategies like household quarantine and the use of anti viral drugs isimpossible for them due to lack of symptoms. However, the result of vaccination is notinfluenced by adding asymptomatic cases. To calculate individual R for the models withboth symptomatic and asymptomatic cases, two (cid:101) p ij values for both the symptomatic andasymptomatic cases have to be calculated using Equation 3. Then, R ,s ( i ) and R ,a ( i )are calculated using the corresponding (cid:101) p ij values. The final R ( i ) values are obtained bytaking the weighted average of R ,s ( i ) and R ,a ( i ).In the study of Longini et al. (2004), a population of 2000 persons in four identicalneighbourhoods is considered. Each individual mixes with people in community, neigh-18 ndividual R0s F r equen cy Figure 1: The frequency of Individual R s without household quarantineborhood, family and play groups. Family sizes differ between one and seven. We havea similar model in the study of Longini et al. (2004) but we also added a mixing groupwork for adults. We constitute a population matrix each row of which includes the ID ofcommunity, neighborhood, family, school-work and the age group of an individual similarto the rows of Table 7. So the number of rows of that population matrix is 2000. The IndividualID FamilyID Size ofFamily NeighborhoodID CommunityID AgeGroup School-WorkID1 10 1 100 1 6 90012 11 2 100 1 6 90013 11 2 100 1 3 3001
Table 7: Population matrix for a model with overlapping mixing groupsdetails of the R code for generating such a population matrix based on census data isavailable from the authors. A major practical problem for this model is the calibrationof the probability of infection within each mixing group. We consider the same infectionprobabilities as in the study of Longini et al. (2004). As disease duration, 3, 4, 5 and6 days with probabilities 0.3, 0.4, 0.2, and 0.1 is again assumed. In the study, we as-sume that an infected person is symptomatic with probability 0.67 and an asymptomaticinfection is only half as infectious as a symptomatic infection (Longini et al. , 2004).In Figure 1, we present the histogram of all individual R values of the population.The figure indicates that only a small number of individuals in the population have19ndividual R values smaller than 1 while most of the population has individual R values between 1 and 2. Moreover, Longini et al. (2004) estimate R as the average of allsecondary cases that the randomly selected initial infective person would infect over allmixing groups he belongs to. They empirically calculate R and find it as 1.7 with a rangeof secondary cases from zero to 17. We compute R ( i ) values with Equation 4 where ˜ p ij is computed by considering the the probability of infection within each mixing group andthe population matrix. Moreover, the average of R ( i ) values for i = 1 , , ... R suggested in the study of Longini et al. (2004) without implementingsimulation. Furthermore, there are many possible R ( i ) values giving the same averageand the importance of R ( i ) values increase as the heterogeneity level increases.In studies with overlapping mixing groups, we found only the suggestion of randomvaccination and of random vaccination of children (Germann et al. , 2006; Longini et al. ,2004). In these studies, the results of vaccination are analysed estimating attack rates viasimulation. However, we suggest to assess the intervention strategies without simulationalso for individual based models. Since we can compute individual R value for eachmember of the population exactly, see the histogram in Figure 1, we can compare thefrequency histograms of individual R values after different intervention strategies areimplemented. Moreover, simulation is not needed while vaccinating individuals based ontheir individual R values by implementing Algorithm 1. In this section, we apply tosimulation only for random vaccination in order to compare intelligent vaccination strat-egy with random vaccination where the vaccinated individual is selected randomly. Wecompare the performance of intelligent vaccination and random vaccination by consider-ing 30%, 50% and 80% of the population vaccinated respectively. Moreover, we recordthe maximum individual R value for the intelligent vaccination. However, we record theminimum, median and maximum of maximum individual R values for random vaccina-tion because each run yields a different maximum individual R value. In Table 8, wepresent the results.Table 8 indicates that 50% random vaccination of the population cannot reduce maxi-mum individual R below 1 in 1000 repetitions while 50% vaccination based on individual R values of the population reduces maximum individual R much lower than 1. Thus,similar to the model with multiple cities and the model with a population partitioned20 accinationPercentage Minimum of MaxIndividual R s in1000 repetitions Median of MaxIndividual R s in1000 repetitions Maximum of MaxIndividual R s in1000 repetitions Max R afterIndividual R BasedVaccination0 4.27 4.27 4.27 4.2730 2.40 3.08 3.73 1.3250 1.70 2.23 2.98 0.9180 0.60 0.96 1.54 0.48
Table 8: Maximum Individual R Values after Random Vaccination and VaccinationBased on Individual R without Household Quarantine Individual R0s F r equen cy Figure 2: The frequency of individual R s after vaccination without household quarantineinto households, we take the advantages of our R formula. Furthermore, the minimumrequired number of vaccinated individuals to guarantee herd immunity is 869. In Figure2, we present individual R values of the unvaccinated population after vaccination of869 individuals based on their individual R s.Finally, we also check how individual R values change if 80% of the symptomaticcases stay home after their first day of infection. Figure 3 shows that even if a consider-able reduction in individual R values is obtained, herd immunity cannot be guaranteedsince one third of the infectious cases are considered to be asymptomatic and house-hold quarantine cannot be implemented for asymptomatic cases. This is also true for100% compliance rate since the actual compliance rate can be at most 67% that is thepercentage of symptomatic cases.The results certainly depend on disease parameters and population structures but we21 ndividual R0s F r equen cy Figure 3: The frequency of individual R s after 80% household quarantineexpect that recursive individual R based vaccination gives consistently a better perfor-mance. So we can see that the calculation of individual R can be a useful tool to assessthe performance of vaccination strategies and also to develop vaccination strategies forstochastic models with arbitrary heterogeneous contact structures. In this paper, we consider a discrete time stochastic SIR model for non homogeneouspopulations and make three main contributions. Firstly, we introduce individual R and propose a general formula for it that is applicable to all types of heterogeneouspopulations with any size. Our other major contribution is the assessment of interventionstrategies by using the formula for individual R without reverting simulation. Lastly,we define intelligent intervention strategies based on individual R values.As we have studied the notion of R for non homogeneous populations, we introducedindividual R as the expected number of secondary cases produced by a unique giveninitially infected individual. In the literature, R for non homogeneous populations iseither calculated by using Markov chains assuming exponential disease time and smallpopulation size or estimated via simulation. Here, we propose a general formula for exactcalculation of individual R that is applicable to an arbitrary mixing structure and largepopulation size.Furthermore, the evaluation of intervention strategies is of practical importance. The22ffectiveness of these strategies is evaluated by simulation studies comparing the averageattack rates or similar characteristics. However, we show that it is possible to assessthe impact of the intervention strategies by using directly the individual R formula. Weanalyze the effectiveness of different intervention strategies by their ability to decrease themaximum individual R value below one. This method is more accurate than descriptivesimulation results to decide how to make an outbreak impossible. However, it is onlypossible to evaluate strategies that are implemented before infection and immediatelyafter one case is infected by using the individual R values of implementing vaccination,household quarantine or the use of antiviral drugs.Finally, an intelligent vaccination policy is developed based on individual R values.Here, the aim is to obtain the greatest reduction in the maximum individual R value fora fixed number of vaccines. It is observed that the number of required individuals to bevaccinated for herd immunity is much higher for random vaccination than vaccinationbased on individual R . References
Ajelli, M., B. Gon¸calves, D. Balcan, V. Colizza, H. Hu, J. J. Ramasco, S. Merler, andA. Vespignani, 2010, “Comparing large-scale computational approaches to epidemicmodeling: agent-based versus structured metapopulation models”,
BMC infectiousdiseases , Vol. 10(1), pp. 190.Allen, L. J., and A. M. Burgin, 2000, “Comparison of deterministic and stochastic SISand SIR models in discrete time”, In:
Mathematical biosciences , Vol. 163(1), pp. 1–33.Antreya, K. B., 2006, “Branching process”,
Encyclopedia of Environmetrics , Vol. 1.Artalejo, J. R. and M. J. Lopez-Herrero, 2013, “On the exact measure of disease spreadin stochastic epidemic models”,
Bulletin of mathematical biology , Vol. 75(7), pp. 1031–1050.Ball, F., 1986, “A unified approach to the distribution of total size and total area underthe trajectory of infectives in epidemic models”,
Advances in Applied Probability , Vol.18(2), pp. 289–310. 23all, F. G., and O. D. Lyne, 2002, “Optimal vaccination policies for stochastic epidemicsamong a population of households”,
Mathematical Biosciences , Vol. 177, pp. 333–354.Barth´elemy, M., C. Godreche, and J. M. Luck, 2010, “Fluctuation effects in metapopula-tion models: percolation and pandemic threshold”,
Journal of theoretical biology , Vol.267(4), pp. 554–564.Bian, L., 2004, “A conceptual framework for an individual-based spatially explicit epi-demiological model”,
Environment and Planning B: Planning and Design , Vol. 31(3),pp. 381–395.Carley, K. M., D. B. Fridsma, E. Casman, A. Yahja, N. Altman, L. Chen, B. Kaminsky,and D. Nave, 2006, “BioWar: scalable agent-based model of bioattacks”, In:
IEEETransactions on Systems, Man, and Cybernetics-Part A: Systems and Humans , Vol.36(2), pp. 252–265.Carvalho, S. A., S. O. da Silva, and I. da Cunha Charret, 2019, “Mathematical model-ing of dengue epidemic: control methods and vaccination strategies”, In:
Theory inBiosciences , Vol. 138(2), pp. 223–239.Colizza, V., and A. Vespignani, 2008, “Epidemic modeling in metapopulation systemswith heterogeneous coupling pattern: Theory and simulations”, In:
Journal of theo-retical biology , Vol. 251(3), pp. 450–467.Craft, M. E., H. L. Beyer, and D. T. Haydon, 2013, “Estimating the probability of amajor outbreak from the timing of early cases: an indeterminate problem?”, In:
PLoSone , Vol. 8(3), pp. e57878.Diekmann, O., J. A. P. Heesterbeek, and J. A. Metz, 1990, “On the definition and thecomputation of the basic reproduction ratio R 0 in models for infectious diseases inheterogeneous populations”,
Journal of mathematical biology , Vol. 28(4), pp. 365–382.Dimitrov, N. B., and L. A. Meyers, 2010, “Mathematical approaches to infectious diseaseprediction and control”, In:
Risk and optimization in an uncertain world , INFORMS.24conomou, A., A. G´omez-Corral, and M. L´opez-Garc´ıa, 2015, “A stochastic SIS epi-demic model with heterogeneous contacts”,
Physica A: Statistical Mechanics and itsApplications , Vol. 421, pp. 78–97.Germann, T. C., K. Kadau, I. M. Longini, and C. Macken, 2006, “Mitigation strategiesfor pandemic influenza in the United State”,
Proceedings of the National Academy ofSciences , Vol. 103(15), pp. 5935–5940.Hernandez-Suarez, C. M., 2002, “A Markov chain approach to calculate R0 in stochasticepidemic models”,
Journal of theoretical biology , Vol. 215(1), pp. 83–93.Inaba, H., 2012, “On a new perspective of the basic reproduction number in heteroge-neous environments”,
Journal of mathematical biology , Vol. 65(2), pp. 309–348.Keegan, L. T., and J. Dushoff, 2016, “Estimating finite-population reproductive numbersin heterogeneous populations”,
Journal of theoretical biology , Vol. 397, pp. 1–12.Kumar, A. and K. Goel, 2020, “A deterministic time-delayed SIR epidemic model: math-ematical modeling and analysis”,
Theory in Biosciences , Vol. 139(1), pp. 67–76.Laskowski, M., B. C. Demianyk, J. Witt, S. N. Mukhi, M. R. Friesen, and R. D. R.McLeod, 2011,
IEEE Transactions on Information Technology in Biomedicine , Vol.15(6), pp. 877–889.Levins, R., 1968,
Evolution in changing environments: some theoretical explorations ,Princeton University Press, New Jersey.Lipsitch, M., T. Cohen, B. Cooper, J. M. Robins, S. Ma, L. James, G. Gopalakrishna, S.K. Chew, C. C. Tan, M. H. Samore, D. Fisman, and M. Murray, 2003, “Transmissiondynamics and control of severe acute respiratory syndrome”,
Science , Vol. 300(5627),pp. 1966–1970.Longini, I. M., M. E. Halloran, A. Nizam, and Y. Yang, 2004, “Containing pandemicinfluenza with antiviral agents”,
American journal of epidemiology , Vol. 159(7), pp.623–633. 25ongini, I. M., A. Nizam, S. Xu, K. Ungchusak, W. Hanshaoworakul, D. A. T. Cummings,and E. Halloran, 2005, “Containing pandemic influenza at the source”,
Science , Vol.309(5737), pp. 1083–1087.L´opez-Garc´ıa, M., 2016, “Stochastic descriptors in an SIR epidemic model for heteroge-neous individuals in small networks”,
Mathematical biosciences , Vol. 271, pp. 42–61.Lund H., L. Lizana and I. Simonsen, 2013, “Effects of city-size heterogeneity on epidemicspreading in a metapopulation: a reaction-diffusion approach”,
Journal of statisticalphysics , Vol. 151(1-2), pp. 367–382.Meyers, L., 2007, “Contact network epidemiology: Bond percolation applied to infectiousdisease prediction and control”,
Bulletin of the American Mathematical Society , Vol.44(1), pp. 63–86.Ross, J. V., 2011, “Invasion of infectious diseases in finite homogeneous populations”,
Journal of theoretical biology , Vol. 289, pp. 83–89.Sharomi, O., and T. Malik, 2017, “Optimal control in epidemiology”,
Annals of Opera-tions Research , Vol. 251(1-2), pp. 55–71.Trapman, P., 2007, “On analytical approaches to epidemics on networks”,
Theoreticalpopulation biology , Vol. 71(2), pp. 160–173.Watson, R., 1980, “A useful random time-scale transformation for the standard epidemicmodel”,
Journal of Applied Probability , Vol. 17(2), pp. 324–332.Wu, J. T., S. Riley, C. Fraser, and G. M. Leung, 2006, “Reducing the impact of the nextinfluenza pandemic using household-based public health interventions”,
PloS medicine ,Vol. 3(9), pp. e361.Yang, Y., P. Atkinson, and D. Ettema, 2008, “Individual space–time activity-based mod-elling of infectious disease transmission within a city”,