The impact of household structure on disease-induced herd immunity

The disease-induced herd immunity level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_D$$\end{document}hD is the fraction of the population that must be infected by an epidemic to ensure that a new epidemic among the remaining susceptible population is not supercritical. For a homogeneously mixing population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_D$$\end{document}hD equals the classical herd immunity level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_C$$\end{document}hC, which is the fraction of the population that must be vaccinated in advance of an epidemic so that the epidemic is not supercritical. For most forms of heterogeneous mixing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_D<h_C$$\end{document}hD<hC, sometimes dramatically so. For an SEIR (susceptible \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document}→ exposed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document}→ infective \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document}→ recovered) model of an epidemic among a population that is partitioned into households, in which individuals mix uniformly within households and, in addition, uniformly at a much lower rate in the population at large, we show that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_D>h_C$$\end{document}hD>hC unless variability in the household size distribution is sufficiently large. Thus, introducing household structure into a model typically has the opposite effect on disease-induced herd immunity than most other forms of population heterogeneity. We reach this conclusion by considering an approximation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{h}}_D$$\end{document}h~D of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_D$$\end{document}hD, supported by numerical studies using real-world household size distributions. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=2, 3$$\end{document}n=2,3, we prove that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{h}}_D>h_C$$\end{document}h~D>hC when all households have size n, and conjecture that this inequality holds for any common household size n. We prove results comparing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{h}}_D$$\end{document}h~D and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_C$$\end{document}hC for epidemics which are highly infectious within households, and also for epidemics which are weakly infectious within households.


Introduction
During the ongoing COVID-19 pandemic there has been considerable discussion of herd immunity.For a very wide range of epidemic models, specifically models for which the basic reproduction number R 0 is given by the maximal eigenvalue of a nextgeneration matrix, if R 0 is greater than one, then vaccinating a fraction h C = 1 − R −1 0 of the population, chosen uniformly at random, with a perfect vaccine (i.e. one that necessarily renders its recipient immune to the disease) in advance of an outbreak reduces the reproduction number to one and thus prevents a large outbreak (see, for example, Diekmann et al. (2013), page 199).The quantity h C is the classical (or vaccination-induced) herd immunity level.For a disease in which infection confers immunity to subsequent infection, herd immunity can also be attained by letting an epidemic run its natural course, possibly with some restrictions in place, for example, lockdown or other non-pharmaceutical interventions.The disease-induced herd immunity level h D is the fraction of the population that needs to be infected before the effective basic reproduction number (i.e.R 0 for an epidemic among the remaining susceptible individuals) is reduced to one.For definiteness, we define h D assuming no restrictions are in place and the epidemic simply runs its natural course.
For an epidemic among a homogeneously mixing population, the classical and disease-induced herd immunity levels are equal.However, that typically is not the case for epidemics among heterogeneous populations.For example, Britton et al. (2020) showed that in a model for COVID-19 in which the population was structured by age and activity level, when R 0 = 2.5, the disease-induced herd immunity level h D = 0.43, which is substantially lower than h C = 0.6, and Gomes et al. (2022) obtained even lower values for h D in a model where individuals varied in susceptibility.These observations have a simple intuitive explanation.For example, in the model with varying susceptibility, individuals with higher susceptibility are more likely to be infected early in the epidemic and consequently the average susceptibility of the remaining susceptibles decreases as the epidemic progresses leading to h D < h C .It seems likely that similar arguments hold for many other forms of heterogeneities, with the general conclusion that introducing heterogeneity into a model has the effect of reducing the disease-induced herd immunity level h D .
An important population structure for epidemics among human populations, which can have a significant impact on disease dynamics and the performance of vaccination strategies, is that induced by households.The aim of this paper is to investigate the impact of household structure on the disease-induced herd immunity level.We use an extension of the SIR (susceptible → infective → recovered) model introduced by Ball et al. (1997) to include an exposed (latent) period.In this model, infective individuals make two types of infectious contacts: local contacts with individuals chosen uniformly at random from their household and, at a much lower rate, global contacts with individuals chosen uniformly at random from the whole population.In sharp contrast to most other forms of heterogeneity, we find that the effect of household structure is generally to increase h D .
For most models it is difficult to calculate h D as it requires knowledge of the trajectory of the epidemic, which typically is not available in closed form.Moreover, in a stochastic model the disease-induced herd immunity level is in fact a random variable, as it depends on the realised trajectory, which converges to h D as the population size converges to infinity.The following approximation to h D , which we adapt to our model, is used in Britton et al. (2020) in the context of a multitype SIR epidemic model.A new model is considered in which all transmission rates are multiplied by a factor κ < 1 and its (limiting deterministic) final outcome is determined.(Note that this factor is denoted by α in Britton et al. (2020).)Let κ be the value of κ so that the effective R 0 among the remaining susceptibles is one.Then the fraction of the population ultimately infected by the epidemic with κ = κ gives an approximation to h D , which we denote by h D .Note that h D is not affected by the introduction of a latent period into the model, as the distribution of the final size of a stochastic SEIR (susceptible → exposed → infective → recovered) model is invariant to very general assumptions concerning the latent period.We adopt a similar approach to obtain an approximation h D to h D for the above households model, except that only global transmission rates are multiplied by κ, with local transmission rates unchanged.
The above definition of h D assumes that no restrictions are in place.Let ĥ D be a generic notation for the disease-induced herd immunity level under restrictions.In practice, ĥ D may depend on the precise pattern of restrictions imposed prior to herd immunity being reached (see, for example, Di Lauro et al. (2021)).A commonly-made assumption in modelling restrictions is that at time t ≥ 0 all transmission rates are multiplied by a factor κ(t).Under this assumption, Britton et al. (2021) show that ĥ D is independent of {κ(t) : t ≥ 0} if mixing is separable.Moreover, for the examples in Britton et al. (2020Britton et al. ( , 2021)), numerical studies showed that the precise timings of restrictions had minimal, if any, effect on ĥ D .The situation is more subtle if restrictions are not applied uniformly, where for some models ĥ D can be highly dependent on the pattern of restrictions (Di Lauro et al. 2021).However, numerical studies indicate that is not the case for the present households model, with restrictions affecting only global transmission rates.Note that h D = ĥ D when such restrictions with factor κ are applied throughout the epidemic.Numerical studies suggest that, under many restrictions, h D is a better approximation than h D to ĥ D .
The usual definition of R 0 as the maximal eigenvalue of a next-generation matrix, or in non-mathematical terms as the mean number of infectious contacts made by a typical infective in an otherwise susceptible population, does not hold for the present households SEIR model, since even in the early stages of an epidemic there are likely to be repeat local contacts within a household.Pellis et al. (2012) give an alternative definition of R 0 via a linear approximation of the early phase of an epidemic in terms of generations of infections, which coincides with the usual definition when it is applicable but can also be extended to models with small mixing groups such as the households SEIR model.Calculation of R 0 for the households SEIR model is quite complex.A simpler to calculate reproduction number for the households SEIR model is R * (see Ball et al. 1997), which is based on the proliferation of infected households (rather than infected individuals).Precise definitions of R 0 and R * and discussion of their properties are given in Sect.2.2.The reproduction number R * is useful for determining herd immunity levels, owing to its ease of calculation.However, it is not comparable between different household structures, unlike R 0 , because it is a households-based rather than an individual-based reproduction number.
Before describing the main results of the paper, we need some more notation.Let λ L denote the individual-to-individual local infection rate and λ G denote the overall rate that an infective makes global contacts.Further, let H be a random variable describing the size of a household chosen uniformly at random and H , the size-biased version of H , be a random variable describing the size of the household to which an individual chosen uniformly at random from the population belongs (see Sect. 2.1). Let H denote the ith factorial moment of H and be the ith factorial moment of a geometric distribution with success probability μ −1 H .The complexities of the households model render analytical results comparing h D and h C hard to obtain in general.First, we consider epidemics which are highly locally infectious, in that if one individual in the household becomes infected then the whole household becomes infected.This assumption, which was introduced in Becker and Dietz (1995), is obtained by letting λ L = ∞ in our model.Under this assumption, the following are our main results.
geometric distribution, so H follows a logarithmic distribution.• Theorem 4.4.Suppose that the epidemic is only just above threshold (i.e.R * is only just above 1) and l * = inf k≥2 {k : • Theorem 4.6 gives an ordering of h D and h C for epidemics which are both highly locally and highly globally infectious.The result is not given explicitly here as it requires appreciable further notation.• Theorem 4.7.Suppose that n > 1 and P( H = n) = p = 1 − P( H = 1), where 0 < p < 1, so a fraction p of individuals reside in households of size n and the remainder in households of size 1.
).An expression for λ * G (n, p) involving the root of an algebraic equation is given in Theorem 4.7.If p is close to 1 (i.e. the households are nearly all of size n) then λ G , and thus R 0 , must be exceedingly large in order to obtain h D < h C .
Analysis is also possible for epidemics that are weakly locally infectious, i.e. when λ L is close to 0. We assume without loss of generality that the infectious period T I has mean 1, and write h D (λ L ) and h C (λ L ) to show explicitly the dependence of these herd immunity levels on λ L .Note that the model reduces to a standard homogeneously mixing SEIR epidemic when λ L = 0, for which for all sufficiently small λ L > 0. • Corollary 4.10.Suppose all households are the same size n > 1.Then, for all sufficiently small λ L > 0, we have The final theorem concerns the case when 0 < λ L < ∞ and all households have the same size n.
• Theorem 4.11.For a common household size n = 2 or n = 3, and for any λ G and We conjecture, supported by numerical studies, that Theorem 4.11 holds for all n > 1.
The remainder of the paper is structured as follows.In Sect.2, we define the stochastic SEIR households model underlying our analysis, describe briefly its threshold behaviour, calculation of the reproduction numbers, R * and R 0 , and of the final outcome in the event of an epidemic taking off.We also present a deterministic model which approximates epidemics that take off.In Sect.3, we describe calculation of the vaccine-induced herd immunity level h C , discuss the definition of the disease-induced herd immunity level h D and describe in detail its approximation h D .Theorems concerning comparison of h D and h C are given in Sect.4, with some of the longer proofs being deferred to an appendix.Numerical comparisons of herd immunity levels are given in Sect.5, including illustration of theorems and study of herd immunity levels for real-world household size distributions.In Sect.6, we give some concluding comments and discuss possible directions for future research.

Model definition
We consider an SEIR (susceptible → exposed → infective → recovered) model for an epidemic among a closed and finite population separated into households.This is similar to the model in Ball et al. (1997) with an extra (exposed) state.The household structure is given as follows.We suppose, for n = 1, 2, . . ., there are m n households of size n.There are m = ∞ n=1 m n households (with m < ∞), and the total population size is N = ∞ n=1 nm n .The epidemic begins at time t = 0 with one initial infective (chosen uniformly at random from the population) and with all other members of the population susceptible.
When a given susceptible is contacted by an infective, they become exposed (latent infective) for a time that is distributed according to a non-negative random variable T E , with an arbitrary but specified distribution which is almost surely finite.When their exposed period ends, an individual becomes infectious for a time distributed according to a non-negative random variable T I , with an arbitrary but specified distribution having finite mean.During their infectious period any given infective makes global contacts with any given susceptible according to a Poisson process with rate λ G N .Further, any given infective makes local contacts with any given susceptibles member of their household according to a Poisson process with rate λ L .Once their infectious period ends, an infective recovers and has no further role in the epidemic.When there are no infectives or exposed infectives remaining, the epidemic terminates.Finally, all Poisson processes describing infectious contacts (whether or not either or both of the individuals involved are the same), as well as the random variables for exposed and infectious periods, are assumed to be mutually independent.
Many of the results in the paper are based on approximations which become exact in the limit as the number of households m → ∞ in an appropriate fashion.For m be the fraction of households that have size n.Precise conditions for such asymptotic results are beyond the scope of this paper.We assume that lim m→∞ α To ease the presentation we suppress the dependence on m of parameters such as α (m) n and just use their asymptotic values.

Threshold behaviour
Suppose that the number of households (m) is large.Since the epidemic begins with one initial infective, the probability that an individual contacted globally belongs to a previously infected household is small during the early stages of the epidemic.Thus the early stages of the epidemic can be approximated by a branching process, describing the proliferation of infected households, in which every global contact is with a previously uninfected household (Ball 1996).The offspring mean R * of this branching process, i.e. the expected number of global contacts occurring from a typical contacted household, is a threshold parameter for the households model.Standard branching process theory implies that, in the limit as m → ∞, the epidemic takes off with strictly positive probability if and only if R * > 1.In the event the epidemic takes off, a non-negligible fraction (a large number of households) of the population becomes infected.
The derivation of R * is as follows.For n = 1, 2, . . ., let αn = nm n N be the probability an individual chosen uniformly at random from the population resides in a household of size n.Consider a globally contacted individual in an otherwise fully susceptible household of size n.This individual begins a local outbreak within their household with dynamics determined by local infection since, in the branching process, global contacts are with previously fully susceptible households.Let μ n (λ L ) denote the mean size, including the initial infective, of a single-household epidemic with n members and local infection rate λ L with only the initial infective infected globally.Global contacts from a given individual occur at rate λ G , and such an individual has mean infectious period E[T I ].Wald's identity for epidemics (see Ball 1986, Theorem 2.1) then gives the mean number of global contacts from a given contacted household of size n to be μ n (λ L )λ G E[T I ].By conditioning on the size of a household of a contacted individual, we have that (2.1) In Ball (1986) it is shown that where φ(θ) = E[e −θ T I ] and β k (λ L ) (k = 1, 2, . . . ) are defined recursively by As explained in Sect. 1, a drawback of R * is that it is not comparable between models with different household structures.An alternative threshold parameter, which does not suffer from that defect, is the basic reproduction number R 0 .As also explained in Sect. 1, the usual definition of R 0 does not hold for households models.Instead, R 0 can be defined by considering generations of infectives, via a directed graph associated with an epidemic (Pellis et al. 2012 andBall et al. 2016).Such a graph is constructed by having population members as the vertices.If, during their infectious period, individual x would contact x , a directed edge is drawn from x to x .The initial infective is the only member of generation 0. The generation of a given individual x is the shortest path length from the initial infective to x.Note that this may not coincide with real-time generations of infectives.Then R 0 is defined by the limit, as the population size goes to infinity, of the asymptotic geometric growth rate of the mean generation size (for a full definition, see Ball et al. (2016), Section 1).For k = 1, 2, . . .and i = 0, 1, . . ., k − 1, let μ (k) i be the mean size of the i th generation for an epidemic in a household of size k with 1 initial infective.(For the present SEIR model, these quantities can be computed using methods described in appendix A of Pellis et al. (2012).)Then R 0 in the households model is the unique positive solution λ of where Ball et al. (2016), Section 2.2 for details.Note that, like R * , the critical value of R 0 is 1.More precisely, R 0 = 1 if and only if R * = 1; R 0 > 1 if and only if R * > 1; and R 0 < 1 if and only if R * < 1.We use R * for calculating or proving results pertaining to herd immunity levels, as it is far simpler to determine than R 0 .We use R 0 when making comparisons between models, owing to its improved interpretability.

Final outcome
Consider an epidemic initiated by one initial infective in a population of size N (m households).Let Z m denote the proportion of individuals infected in the epidemic.Provided R * > 1, as m → ∞, Z m converges in probability to a discrete random variable Z with probability mass function for some 0 < z < 1 defined below.(The mass at zero in the random variable Z corresponds to the branching process in the previous section going extinct.)We define a major outbreak to have occurred if Z m ≈ z and it follows that the sum of the infectious periods of all infected individuals in a major outbreak, S m , is approximately Hence, the probability that a randomly chosen individual avoids global infection during the course of a major outbreak is approximately π = exp(−λ G E[T I ]z), since to avoid global infection there must be no points in a Poisson process of intensity λ G /N run for time Let μn (λ L , π) denote the mean size of a single-household epidemic in a household of size n with local infection rate λ L and Bin(n, 1 − π) initial (globally infected) infectives; using the standard Bin(n, p) notation to denote the binomial distribution.Denote this epidemic model, which is considered in Addy et al. (1991), by Ẽn (λ L , π).
Returning to the households model, suppose that a major outbreak occurs and let Tn denote the total number of individuals infected in a typical household of size n, all of whom are initially susceptible.In the limit as m → ∞, individuals independently avoid global infection with probability π .Thus in a household of size n, Bin(n, 1 − π) will be infected globally and hence E[ Tn ] = μn (λ L , π).Then equation (3.10) of Ball et al. (1997) yields The probability that a given individual in an initially susceptible household of size n is infected during the epidemic is approximately μn (λ L , π)/n.Conditioning on the household size of a randomly chosen individual then establishes (2.4) Note that for large m the effect of the atypical behaviour of the household containing the initial infective becomes negligible and disappears in the limit as m → ∞.Thus, since π = exp(−λ G E[T I ]z), (2.4) admits an implicit equation for z, whereby z = 0 is always a solution and a second solution z ∈ (0, 1) exists if and only if R * > 1.
A similar argument can be used to establish, in the event of a major outbreak, the proportions (P n,v ) (n = 1, 2, . . .; v = 0, 1, . . ., n) of households of size n with v members ultimately infected.For n = 1, 2, . .., (P n,v ) satisfies the system of equations (see Addy et al. (1991), equation 4) (2.5) The above arguments are made rigorous in Ball et al. (1997), Section 4.2 and hold with or without the inclusion of a latent period, see Ball et al. (1997), Section 3.1.

Deterministic model
In House and Keeling (2008), Section 2, a system of ordinary differential equations (ODEs) is derived for the evolution over time of the SIR epidemic model with households, assuming that T I ∼ Exp(γ ), i.e. the infectious period distribution is exponential with mean γ −1 .House and Keeling's ODEs represent the deterministic limit of the stochastic process defined in Sect.2.1 as m → ∞, under the assumptions that all households are the same size and there is no latent period.We extend the system of ODEs to allow for variable household size and a latent period T E ∼ Exp(δ).Consider the model in Sect.2.1 with maximum household size n max .Let where Z + = {0, 1, . . .}, and for t ≥ 0 and (s, e, i, r ) ∈ H (n max ) , denote by s,e,i,r (t) the number of households with s susceptible, e exposed, i infectious and r recovered members at time t.For n = 1, 2, . . ., n max , let (2.7) be the set of states in which there is at least one infective or exposed individual.For (s, e, i, r ) ∈ H (n max ) , we assume that 1 m H (m) s,e,i,r (0) → h s,e,i,r (0) as m → ∞, where (s,e,i,r )∈H (nmax) + h s,e,i,r (0) > 0, so a strictly positive fraction of the population is initially either exposed or infective in the limit as m → ∞.Then, under the Markovian assumption, we have that 1 m H s,e,i,r (t) converges in probability to a deterministic process h s,e,i,r (t) as m → ∞; see Ethier and Kurtz (1986), Theorem 11.2.1.Clearly, for n = 1, 2, . . ., n max , we have (s,e,i,r )∈H n h s,e,i,r (t) = α n for all t.Let ī(t) = (s,e,i,r )∈H (nmax) ih s,e,i,r (t).The deterministic model can be obtained by considering the possible transition rates between household states and yields, for (s, e, i, r where, on the right-hand side, h s ,e ,i ,r (t

Herd immunity in SEIR households model
In this section we outline the various versions of herd immunity that we consider.We start by recapping vaccine-induced herd immunity; we then describe and explore the details of the disease-induced herd immunity level h D and its approximation h D which are outlined in Sect. 1.In this section, and throughout the remainder of the manuscript, we assume that R * > 1, since if this is not the case then herd immunity is already achieved.

Vaccine-induced herd immunity level h C
Suppose some members of the population are vaccinated before an epidemic occurs.
Assume that such a vaccine is given to each member of the population independently with probability c, and is perfect, so that vaccinated individuals are completely immune to infection.Then, for v = 0, 1, . . ., n, a given household of size n has v members vaccinated according to the (binomial) We obtain a postvaccination threshold parameter RU (c) by considering a branching process of potential global contacts.If a potential global contact is with a susceptible individual then it triggers a local epidemic; if the contact is with a vaccinated individual then nothing happens.
The mean number of potential global contacts emanating from a single-household epidemic for a household in state (n, v) that is contacted globally, with the initial infective chosen uniformly at random from members of the household (μ n,v , say) is given by for n = 1, 2, . . .and v = 0, 1, . . ., n.This is because a vaccinated member being contacted leads to no global contacts at all, and an unvaccinated member being contacted initiates a single-household epidemic amongst the n −v non-immune members.Such an unvaccinated individual is contacted with probability n−v n .Conditioning on the vaccination status of an individual's household, as well as the individual's household size, and using the same argument as for equation (2.1) yields a post-vaccination threshold parameter with μ 0 (λ L ) = 0.The function RU (c) is continuous and strictly decreasing, with The function g is strictly increasing and X is stochastically decreasing in c.Thus, f n (c) is strictly decreasing in c, whence so is RU (c).)This implies there is a critical value, h C say, such that RU (h C ) = 1 and a major outbreak can be avoided.The quantity h C is the vaccine-induced herd immunity level.Note that by ensuring RU (c) ≤ 1, the whole population is considered protected as a major outbreak is no longer possible.This argument regarding uniform vaccination is considered (among other vaccination strategies) in Ball and Lyne (2006).
As noted at the start of Sect. 1, for epidemic models in which R 0 is given by the maximal eigenvalue of a next generation matrix, 0 for all n max ; see Remark 2 following the aforementioned Theorem 1.If λ L = 0, the model reduces to a standard homogeneously mixing SEIR epidemic and h C = 1 − R −1 0 .

Limiting disease-induced herd immunity level h D
An alternative method of achieving herd immunity in a population arises from the spread of a first wave of infection, in which infected members from the first wave are considered immunised thereafter.
Consider the SIR version of the households model described in Sect.2.1, with T I ∼ Exp(γ ) and assume this epidemic takes off.As the epidemic progresses some members of the population are infected, lowering the overall susceptibility of the population.Suppose that the first epidemic is stopped (i.e.all infectious spread, including that within households, is stopped) at time t > 0. Consider a second epidemic initiated at time t with one initial infective and all those individuals infected by time t in the first epidemic immune to infection in the second.Recalling that m is the number of households in the population, let R (m) V (t) be the threshold parameter (R * ) for this second epidemic, which is a random variable owing to its dependence on the trajectory of the first epidemic.If R (m) V (t) ≤ 1, then the second epidemic is not supercritical and a major outbreak cannot occur.Disease-induced herd immunity is achieved when the trajectory of R (m) 123 denote the fraction of the population that is still susceptible at time t in the first epidemic by S (m) (t).Then the disease-induced herd immunity level and is also a random variable determined by the trajectory of the first epidemic.We conjecture that H (m) , where h D is a constant, and in the presence of a latent period Moreover, we conjecture that, under suitable conditions, these convergence results hold also when T I and T E follow non-exponential distributions.
To compute h L D when T I ∼ Exp(γ ) and T E ∼ Exp(δ), we use the deterministic model in Sect.2.4.Recall the definitions of H (n max ) and H and suppose that h(0) = .Then the deterministic model given by (2.8) can be used in the obvious fashion to define a disease-induced herd immunity level Other than in Appendix A, when computing h L D , we assume that initially a fraction = 10 −5 of households are in state (n max − 1, 0, 1, 0), with all other households being fully susceptible.An equivalent assumption is made when computing h D .Note that, if the largest household is of size n, the system in (2.8) contains exact order n 4 /24 equations (or n 3 /6 in the SIR case) and becomes computationally expensive to solve when n is large.
The proofs of the above conjectures require extending the theory of Barbour and Reinert (2013) to the households model and are beyond the scope of this paper.Some brief comments, together with numerical evidence in support of these conjectures are given in Appendix A.

Approximate disease-induced herd immunity level hD
Calculation of h D requires deterministic limiting equations, which are not tractable in general; in the Markovian setting, for example, h D can be found numerically.For other infectious period distributions, such a calculation is generally not available; we consider an approximation ( h D ) to h D .We adapt the approach taken in Britton et al. (2020) to the households setting, using the final outcome of an epidemic with reduced global infection rate to approximate the state of the population at the time herd immunity is achieved.The method of approximation is as follows: Let κ ∈ R −1 * , 1 .Run to its conclusion an epidemic in which the global infection parameter λ G is replaced by κλ G and all other parameters are unchanged, i.e. an epidemic with prevention measure applied to global infection only.Then expose the population to a second epidemic with κ = 1, with members infected in the first epidemic now immune.The threshold parameter R * for this second epidemic is a function of κ, which we denote by RDI (κ).Determine κ, the smallest value of κ such that RDI (κ) ≤ 1 by solving RDI ( κ) = 1.Then h D is the fraction of the population infected in the first epidemic with κ = κ.In summary, we adjust the global infection rate in the first epidemic to force criticality in the second, and then consider the final size of the first Fig. 1 Distribution of the number of susceptibles in a typical household when herd immunity is achieved under h D and h D when all households have size 4.For each p L , λ G is chosen so that R 0 = 2 epidemic.Note that this method relies on final size results, which are more amenable to study than time-dependent results.

Accuracy of the approximation of h D by hD
The disparity between h D and its approximation h D depends on the distribution of susceptibles inside households when herd immunity is achieved.When λ L = 0, there is no within-household spread; this distribution is the same under h D as under h D .The same conclusion holds when λ L = ∞, since all single-household epidemics end immediately (everyone in a contacted household becomes infected as soon as that household is contacted).Thus, when λ L = 0 or λ L = ∞, we have h D = h D .For epidemics with 0 < λ L < ∞, local epidemics are run to termination under h D , but not under h D .This, coupled with a lower global infection rate under h D , leads to h D being more strongly governed by within-household spread.A consequence of this is the distribution of susceptibles among households being more clumped under h D than h D , leading to a larger proportion of households with no susceptible individuals than under h D , which often results in both h D > h D and the approximation becoming worse as n max increases (cf. the discussion following Theorem 4.1 in Sect.4.2.2).The above-mentioned clumping is illustrated in Fig. 1, which considers the case of a common household size n = 4, with p L varying in [0, 1] and λ G being chosen so that R 0 = 2.When p L = 0, there is no within-household spread and the distribution of the number of susceptibles in a typical household when herd immunity is reached is Bin(4, 0.5) under both h D and h D .Note the agreement in the two distributions of susceptibles when p L = 1.For p L ∈ (0, 1), the distribution has greater mass at the extremes 0 and 4 under h D than h D .We suspect that h D > h D holds generally for a common household size n > 1 (with 0 < λ L < ∞) and give numerical examples supporting this claim in Sect.5.1.It is possible (but atypical) for h D < h D to occur, and the difference is small in the cases we have met.An example is a household structure comprised of households of size 1 and n > 1 only; see Fig. 7 in Sect.5.1.
The accuracy of the approximation of h D by h D is explored in Fig. 2, which shows heat maps of the percentage error 100 h D − h D /h D as a function of ( p L , R 0 ) for common household sizes n = 2, 3, 4 and 5. Observe that the percentage errors are all small, increase with n and are greatest for intermediate values of p L .The maximum percentage error as ( p L , R 0 ) varies over [0, 1] × (1, 25] for each choice of n and for some other household size distributions are given in Table 1.Note that the value of p L at which this maximum is attained tends to decrease with mean household size μ H = ∞ n=1 nα n .This may be a consequence of the fact that for fixed p L the fraction infected by a single-household epidemic increases with household size.The maximum percentage errors are small, except for countries with large mean household sizes.Moreover, these are maximum errors and even if they are not small, the error is  The first four rows correspond to a common household of size n.The fifth row corresponds to α1 = α4 = 0.5 and the remaining rows correspond to real-world household size distributions (see Sect. 5.2) small for many choices of parameter values, as illustrated by the n = 5 heat map in Fig. 2. Thus, h D is generally a very good approximation of h D .

Impact of restrictions
Note that the parameter κ corresponds to restrictions being placed on the population which affect only the global infection rate; the severity of such restrictions increases as κ decreases.Since κ is chosen such that the second epidemic is at criticality, κ corresponds to the most severe restrictions that can be placed for the whole duration of the first epidemic, such that the second epidemic is not supercritical; more severe restrictions will leave the second epidemic supercritical, so herd immunity will not be achieved.Recall that ĥ D denotes the disease-induced herd immunity level when restrictions are in place.Numerical investigations suggest that as κ increases from κ to 1, ĥ D decreases; see Fig. 3.In this example, which uses the UK and Morocco household size distributions (see Sect. 5.2), such restrictions have only a small effect on the herd immunity level ĥ D , with the effect being larger for Morocco.Moreover, when λ L is fixed, we observe that κ decreases as R 0 increases.Repeating these calculations for other values of λ L (not shown) reveals similar patterns and suggests that the effect is largest when λ L is around 0.5.All of these observations regarding the UK and Morocco household size distribution are also seen with larger and more variable household size distributions: the effect of κ on ĥ D remains small but becomes slightly larger if the household size distribution is more variable, κ decreases with increasing R 0 and the effect of varying κ seems greatest when λ L ≈ 0.5.

Outline
This section presents results concerning orderings of h D and h C .Since both of these quantities depend only on final outcome properties of the epidemic, they are invariant to the distribution of the latent period and we therefore take T E = 0, corresponding to the SIR setting, in this section.The problem of solving for h D is not analytically tractable when 0 < λ L < ∞.We hence begin with the highly locally infectious case where λ L = ∞, considered in Becker and Dietz (1995), for which a framework for comparison of h D and h C is established in Sect.4.2.1 and explicit progress is made.This is then applied to several household size distributions, beginning with all households being the same size, where h D > h C is established (Theorem 4.1).Further, in Sect.4.2.2 we show that for a common household size n, the maximum of h D − h C as a function of λ G occurs when λ G = 4 (1+n)E[T I ] , corresponding to R 0 = 2 (Theorem 4.2).A necessary and sufficient condition for h D = h C in the highly locally infectious case (Theorem 4.3) is derived in Sect.4.2.3, leading to study of h D and h C for epidemics just above criticality (Theorem 4.4 in Sect.4.2.4), as well as highly locally and globally infectious epidemics (Theorem 4.6 in Sect.4.2.5).In Sect.4.2.6 we derive a result for household structures with only households of size 1 and n > 1 (Theorem 4.7).The weakly locally infectious case λ L → 0 is treated in Sect.4.3 and a condition for h D > h C is derived (Corollary 4.9).
We consider the general case 0 < λ L < ∞, both analytically and numerically, when all households are of the same size.In Sect.4.4, we prove that h D > h C for a common household size n = 2 and n = 3 (Theorem 4.11).We conjecture that h D > h C holds for common household size n ≥ 4; supporting evidence for this conjecture is provided in Sect.5.1.

General framework
In the highly locally infectious case (λ L → ∞) explicit analytical progress is possible, as any infected individual will infect their whole household.We therefore have μ n (λ L ) = n for n = 0, 1, . . . .Using (2.1), we find that n=1 n αn is the mean size of the household of an individual chosen uniformly at random from the population.Thus, The inner sum in (4.1) can be evaluated using the second moment of a Bin(n, 1 − c) random variable.Using the definition of μ H , it follows that h C is given by the unique solution in (0, 1) of the quadratic equation yielding As noted at the end of Sect.3.1, h C = 1 − R −1 0 in the present highly locally infectious case.
Turning to the disease-induced herd immunity level h D , consider the first epidemic with global infection rate κλ G , where κ solves RDI ( κ) = 1 as described in Sect.3.2.2.Let z( κ) be the fraction of the population infected by that epidemic and π = exp(− κλ G E[T I ]z( κ)), the probability that any given susceptible avoids global contact during that epidemic.For n = 1, 2, . . .and v = 0, 1, . . ., n, let x n,v (π ) be the proportion of households with n members which have v infected in this epidemic and thus immune to the second epidemic.
Letting R DI (π ) = RDI ( κ), we have In the highly locally infectious case, a given individual escapes infection if and only if their whole household avoids global infection.Thus x n,0 (π ) = π n , x n,n (π ) = 1 − π n , and x n,v (π ) = 0 for v / ∈ {0, n}.Substitution into (4.4)yields where Combined with (4.2), we have a framework to compare h D and h C .Note, however, that the system given by (4.6) and (4.7) does not always allow closed-form calculation of h D .Note that in this subsection dealing with the highly local infectious case, the distribution of T I only enters our results through its mean E[T I ].

Common household size
Suppose that all households are of size n.When n = 1 the model reduces to the standard homogeneously mixing model, so Note that f H (π ) = π n , so (4.6) and (4.7) yield . (4.9) Theorem 4.1 Consider the highly locally infectious case with common household size n > 1.
Proof The claim is established by subtracting equation (4.8) from (4.9), and letting for ease of exposition, so Expressing explicitly the dependence of h D and h C on x, we obtain (4.10) The result then follows by elementary manipulations of (4.10), since R * > 1 implies x ∈ (0, n).
A heuristic justification for Theorem 4.1 is as follows.In disease-induced herd immunity, after the first epidemic, households contain either 0 or n susceptibles, depending on whether that household was infected.Consider a randomly chosen individual contacted globally in the second epidemic.If this individual is immune, this contact contributes no further infection.Otherwise, they begin an epidemic within their household which, in the highly locally infectious case, will infect all non-immune members.
Hence, under disease-induced herd immunity, the potential for within-spread is as high as possible (the rest of the household is susceptible).Thus disease-induced herd immunity corresponds to the worst possible vaccination strategy for a given coverage, resulting in h D > h C .
A further result provides a link between the highly locally infectious case and R 0 for the households model, in which we treat h D − h C as a function of λ G .

Theorem 4.2 Under the same assumptions as Theorem 4.1, h D − h C has a unique maximum as a function of λ G which is attained when λ
Proof We show that h D (x) − h C (x) has a unique stationary point, which must be a maximum since, from (4.10), h D (x) − h C (x) → 0 as x ↓ 0 and x ↑ n, corresponding to R * → ∞ and R * → 1, respectively, and by Theorem 4.1, h D (x) − h C (x) > 0 for x ∈ (0, n).Then we find the value of x which yields the maximum and compute the corresponding R 0 value.
Ignoring the denominator in (4.10), differentiation with respect to x and equating to 0 leads us to solve In the highly locally infectious case, all secondary infections in a household are attributed to the primary case, so μ Setting λ G = λG in (2.3) gives that R 0 is the unique positive root of g n (λ) = 0, where

Necessary and sufficient condition for hD = h C for all G
The framework given in Sect.4.2.1 for the highly locally infectious case enables proof of the following result.For θ ∈ (0, 1], we write H ∼ Geom(θ ) if H has a geometric distribution with probability mass function

Theorem 4.3 In the highly locally infectious case, h
. We begin by assuming and solving for H . Equations (4.2) and (4.6) give .
Further, f H (1) = 1 since f H (π ) is a probability-generating function.This separable ODE can be solved to yield, for 0 ≤ π ≤ 1, which is precisely the probability-generating function of a Geom(μ −1 H ) random variable.This establishes the only if part of the equivalence claim.For the converse, assume that H follows a geometric distribution with parameter μ H . Then (4.11) holds and the logic for the above proof is reversible thereafter, so the result follows.
When H follows a geometric distribution, H follows a logarithmic distribution.Note that real-life household size distributions will have a finite maximum size, so for any realistic household size distribution h D = h C will not hold for all λ G in the highly locally infectious case.We typically observe h D > h C for real-life household size distributions and comment upon this further in Sect.5.1.

Just supercritical epidemics
We can use the framework provided in Sect.4.2.1 to give an ordering of h D and h C for epidemics which are just above threshold, i.e. when R * is just above 1, so π is just below 1.We hence establish an ordering of h D and h C by considering G(π ), given in (4.11), in the neighbourhood of π = 1.Assuming a fraction z(π ) of individuals are infected in the first epidemic gives a threshold parameter R DI (π ) for the second epidemic, as given in (4.4).Vaccinating the same proportion uniformly at random gives a threshold parameter R U Clearly we have G(1) = 0, since f H (1) = 1 and f H (1) = μ H . Let G (k) be the k th derivative of G and define = 1 and it follows that h D > h C for epidemics which are just above threshold.A similar argument shows that Determining the ordering of h D and h C reduces to comparing factorial moments of H to those of a geometric distribution with parameter μ −1 H .For a random variable H define, for i = 1, 2, . . ., the factorial moments μ [i]  H Theorem 4.4 Let H be a given size-biased household size distribution with mean μ H and factorial moments μ [i]  H (i = 0, 1, . . .).Suppose that l * = inf
Similarly to the just supercritical case, the only distribution for H which has G (n) (0) = 0 for all n is the geometric distribution with parameter μ −1 H . Theorem 4.6 then reduces the ordering of h D and h C to iterative comparison of the probability mass function of H with the relevant geometric distribution.

Households of size 1 and n > 1
Theorem 4.1 in Sect.4.2.2 shows that, in the highly locally infectious case, h D > h C for all λ G when the households all have the same size.We now consider the simplest setting when there is variability in household size, i.e. the case where there are only two household sizes, 1 and n > 1.For 0 < p < 1, let p denote the proportion of individuals who belong to a household of size n.Thus αn = p and α1 = 1 − p.We consider how h D − h C varies with p, with a view towards obtaining different orderings of h D and h C as the household structure changes.The following theorem, proved in Appendix B.1, shows that the ordering of h D and h C is less straightforward when household size is variable.
, where πn ( p) is the unique root in (0, 1) of One can solve (4.17We see immediately that R * 0 (n, p) decreases with n and increases with p.One can also show that, for fixed n, lim p↑1 R * 0 (n, p) = ∞.Further, var( H ) = (n−1) 2 p(1− p), so the variability in household size is small when p is close to one.Thus Theorem 4.7 shows that, even with low variability in household size, we can have h D < h C ; however, in this example, R 0 has to be unrealistically large for this to happen.

Weakly locally infectious case
Consider now the case when the extra local infection is small but non-zero, corresponding to λ L → 0, with R * > 1, E[T I ] = 1 and var(T I ) < ∞.Assume also that n max < ∞.
The proof of Theorem 4.8 involves computing the first three terms of the Maclaurin expansion of h D (λ L ) − h C (λ L ) and is given in Appendix B.2, where the assumption that n max < ∞ is also explained.The assumption that E[T I ] = 1 involves no loss of generality (since time can be rescaled appropriately) and is made to simplify the presentation of the proof.The assumption var(T I ) < ∞ is required for the third term in the above-mentioned Maclaurin expansion.Note that when 0).The following corollary is an immediate consequence of Theorem 4.8.
are required in order to give an ordering.Note that Corollaries 4.5 and 4.9 produce contrasting orderings of h D and h . (See Sect.5.2 for a numerical exploration of this.) A result for a common household size n > 1 follows immediately from Corollary 4.9.
Corollary 4.10 Suppose all households are the same size n > 1.For all sufficiently small λ L > 0, we have Proof When the common household size is n > 1, we have var( H ) = 0 and E[ H − 1] > 0. Applying Corollary 4.9 then establishes the claim.

Common household size with 0 < L < ∞
We have shown that, for a common household size n and any λ G such that R * > 1, we have h D > h C when λ L → 0 and λ L = ∞.The following theorem considers λ L ∈ (0, ∞).Theorem 4.11 For a common household size n = 2 or n = 3, and for any λ G and Proof Suppose a fraction of the population is infected by a first epidemic in the households model with the above parameters.This leads to a threshold parameter RDI (z) for the second epidemic.Further, assuming the same proportion are vaccinated uniformly at random gives threshold parameter RU (z).We show that RDI (z) − RU (z) > 0, from which h D > h C is immediate.
Let P D i (i = 0, 1, . . ., n) be the proportion of households with i members immune owing to the first epidemic and let P U i (i = 0, 1, . . ., n) be the analogous quantity with uniformly at random vaccination, both assuming a fraction z of the population is immune.Then we can write and Assuming n = 2 and considering the proportion of susceptibles remaining after vaccination also yields since the probability an individual avoids global infection π is larger than the overall probability it avoids infection 1 − z.We then find and A > 0. The result follows and the claim is established for n = 2.The proof for n = 3 uses a similar (but more involved) argument and is deferred to Appendix B.3.
A proof for n > 3 has not been forthcoming, but we make the following conjecture, which is supported by numerical evidence (Fig. 6) in Sect.5.1.Conjecture 4.12 For any common household size n > 1, and for any λ G and λ L such that R * > 1, we have h D > h C .
In Theorem 4.2 we show that for a common household size in the highly locally infectious case, the difference h D − h C is maximised at R 0 = 2.We show numerically that this does not hold when λ L < ∞.For ease of visualisation we work in terms of the probability that an infectious individual makes local infectious contact with a given individual in their household, maximises the difference (assumed to be positive on the basis of Conjecture 4.12) between the herd immunity levels, then calculate the resulting value of R 0 .
We see that the 'optimal' value of R 0 broadly increases with p L and tends to 2 as p L → 1, consistent with Theorem 4.2.The dip near p L = 0 becomes more pronounced as n increases.Similar behaviour occurs for other choices of infectious period distribution.
When p L = 0 we have h D = h C , so any value of λ G maximises h D − h C .As such, the value of R 0 when λ L = 0 is not well-defined, leading to instability when solving numerically.However, in the general setting with variable household sizes we can proceed analytically using Theorem 4.8.We have, as λ L ↓ 0, where

Numerical comparisons of herd immunity levels
We begin in Sect.5.1 by numerically illustrating some of the results from Sect. 4, then in Sect.5.2 we explore how our findings play out in the context of realistic household size distributions.Throughout this section, we restrict attention to the Markovian case, i.e. we assume that T I ∼ Exp(γ ), and T E ∼ Exp(δ) if a latent period is present.

Illustrative examples
In this subsection we assume that γ = 1 and, where applicable, δ = 1 also.We begin by plotting in Fig. 6 how h D , h L D , h D and h C vary with p L for a common household size n when λ G is chosen so that R 0 = 2 is fixed.
Figure 6 provides an illustration of Theorem 4.11 and support for its extension in Conjecture 4.12; we observe h D > h C throughout Fig. 6.We also observe from Fig. 6 that we have h D > h L D .Analytical comparison of h L D and h D is often not tractable, however we can show that h D > h L D for p L sufficiently close to 1 as follows.When p L = 1, under the SIR model disease-induced herd immunity leaves households either fully susceptible or fully non-susceptible.As noted in the discussion following Theorem 4.1, this corresponds to the worst possible vaccination strategy for a given coverage, implying h D > h L D , since under the SEIR model there may be households in which only the initial case in that household is non-susceptible.
In Fig. 7 we consider the same comparisons as in Fig. 6, but with a variable household size distribution.Specifically, we take α1 = αn = 0.5 for some n > 1; half the individuals in the population are in households of size 1 and the other half are in households of size n.This implies that α 1 = n/(n + 1) and α n = 1/(n + 1).A first observation based on Fig. 7 is that the behaviour near p L = 0 is consistent with the predictions of Sect.4.3.Specifically, since E[ H ] = (n − 1)/2 and var( H ) = (n − 1) 2 /4, Theorem 4.8 predicts that h C > h D near p L = 0 if and only if n ≥ 4.There is contrasting behaviour in terms of the shape of the herd immunity levels as n increases.When n = 2 and n = 3, h C is the smallest of the considered herd immunity levels for all p L .By contrast, when n = 4 or n = 5 there are values of p L such that h C is the largest of the considered herd immunity levels.As n increases in Fig. 7, the approximation h D for h D gets worse, cf.Sect.3.2.3.Introducing a latent period does not necessarily lead to a reduction in the disease-induced herd immunity level; we observe that h L D > h D when n = 4 and n = 5.Note that in Figs. 6 and 7 we have h C ≥ 1 − 1/R 0 , with strict inequality unless all households are of size 3 or less; this follows from Theorem 1 of Ball et al. (2016).

Real-world household size distributions
This section is motivated by the study in Britton et al. (2020), which considers the influence of population heterogeneity on the disease-induced herd immunity level for the COVID-19 pandemic.Britton et al. (2020) uses a Markov SEIR model in a popula-tion that is structured by age and activity level, in which for all individuals the exposed and infectious periods follow Exp(1/3) and Exp(1/4) distributions, respectively, with the unit of time being a day.Thus, the mean exposed and infectious periods are 3 and 4 days, respectively.Using the approximation to h D described in Sect. 1, Britton et al. (2020) find that, when R 0 = 2.5, h D for a homogeneously mixing model and h C are both 60%; for the model with both age and activity structure, h D is reduced to 43.0%.
The aim of the present numerical study is to investigate the effect of household structure on h D , using a model with the above values of δ and γ and a range of real-world household size distributions.In order to achieve that we need a way of calibrating models with different choices of (λ L , λ G ).One possibility is to keep the basic reproduction number R 0 fixed.However, R 0 is not uniquely defined for household models.The definition in Sect.2.2 uses so-called rank generations and a different value for R 0 would typically be obtained if real-time generations were used instead, as for example in Neal and Theparod (2019).In practice, for an emerging epidemic, R 0 is often estimated indirectly, via an estimate of the epidemic's early exponential growth rate r ; see, for example, Wallinga and Lipsitch (2007).For the multitype SEIR model used in Britton et al. (2020), R 0 and r satisfy see Sections 1.3.1 and 1.5 of the supplementary material of Trapman et al. (2016).Note that the relationship (5.1) between R 0 and r is the same for all models in this class of multitype Markov SEIR epidemics and in particular matches that for the homogeneously mixing Markov SEIR model (Trapman et al. 2016).
We adopt the following method of calibrating models with different (λ L , λ G ), based on the early exponential growth rate r .For a given choice of R 0 in Britton et al.'s model, which we denote by R BBT The region enclosed between the solid and dashed black curves in Fig. 8 represents the set of values of (μ H , σ 2 H ) for which Corollary 4.5 and Corollary 4.9 give different orderings of h D and h C .We see that all countries considered have household size distributions in this set; though some are very close to the critical line E[ H − 1] = var( H ) in the weakly locally infectious case.
We now explore the various herd immunity levels in our SIR and SEIR models, using the household size distributions of the UK (Fig. 9) and Morocco (Fig. 10) as exemplars.These countries are chosen because of their quite different household size distributions (cf.Fig. 8).The computation of h L D is omitted for Morocco as its calculation becomes numerically infeasible, since the dimension of the system of ODEs (2.8) becomes too large owing to the high maximum household size.
Considering the UK household size distribution, which has μ H = 3.02 and σ 2 H = 2.26, we see that h D > h C , which is as expected given we have observed h D < h C only in cases where household sizes have very high variability.We also observe less variation in h L D than in the other herd immunity levels.Increasing R BBT 0 leads to the growth rate r being fixed at a higher value, in turn causing higher herd immunity levels.We also observe that h L D and h D are very close for fixed r as p L increases from 0, until around p L = 0.6.
We observe very similar qualitative behaviour for other household size distributions, as shown for the Morocco household size distribution (which has μ H = 5.74 and σ 2 H = 6.12) in Fig. 10.We now explore the quantitative differences between the herd immunity levels in detail for a wider range of countries' household size distributions.Specifically, in Fig. 11 we compare the various herd immunity levels between several countries in the absence of a latent period, with R BBT 0 = 3. (Estimates of R 0 for COVID-19 vary greatly even for the same country, but other choices for R BBT 0 produce qualitatively similar results.) We observe h D > h D in Fig. 11, as well as h D = h D at p L = 0 and p L = 1.For countries with generally smaller household sizes (i.e.Finland, Japan and the UK), h D and h D are very close in value.Countries with a larger value of μ H give larger values for h D and h D but lower values of h C .We generally observe h D > h C ; the exceptions to this are for Morocco, Finland and Japan when p L is close to zero, and even then the difference between h D and h C is very small.We see h C decreases monotonically with p L , whereas h D and h D are not monotone in their dependence on p L .Finally, considering the last plot in Fig. 11, we find that the difference h D − h D is maximised at a smaller value of p L , with a larger maximum difference, when μ H is larger.

Concluding comments
We have presented a general framework for investigating disease-induced herd immunity in epidemic models with household structure.Calculating the disease-induced herd immunity level h D for such models is not straightforward and we have introduced a useful approximation h D , which is more amenable to analysis.In sharp contrast to most forms of heterogeneous mixing, for which h D is less than the vaccine-induced herd immunity level h C , the imposition of household structure generally leads to h D being greater than h C , unless the variability in the household size distribution is suffi- ciently large.This is proved using h D for epidemics which are either highly or weakly locally infectious, and numerical studies support the conjecture that it holds more generally.
It would be worthwhile to consider more fully the impact of restrictions, such as lockdown, on ĥ D , the disease-induced herd immunity level when restrictions are in place.In Sect.3.2.2,we give an example where such restrictions affect only the global infection parameter λ G , for which the impact of the restrictions on ĥ D is minimal; moreover, the approximation of ĥ D by h D improved with increasing restrictions.Similar results were found with other examples.However, in that example restrictions were applied uniformly with time which is unlikely to be the case in practice.Further, in practice restrictions may also affect the local infection rate λ L ; indeed it is not hard to envisage scenarios in which λ L might increase.
Another worthwhile avenue for future research is to consider models which combine household structure with other forms of heterogeneous mixing.This can be done using the multitype households model and a similar approximation to h D for the diseaseinduced herd immunity level can be calculated using results in Ball and Lyne (2001).We are currently investigating this for a model with activity levels, as in Britton et al. (2020), and also household structure.We have assumed throughout that the individual-to-individual local infection rate λ L is independent of household size n.Although this assumption is often made with household models and is usually reasonable for small n, such as in the UK, Sweden and Finland household size distributions, it is less easily justified for countries with large household sizes, such as Iraq, Pakistan and Chad.One would expect λ L to decrease with n and it would be interesting to explore the consequent impact on h D .Note that the results of Sect.4.2 concerning the highly locally infectious case are unaffected but other results may change.
Throughout a large part of this work we have used h D as an approximation to h D .The only models in which we have computed h D are those in which the infectious and latent periods follow exponential distributions.In real-life epidemics, the distributions of these quantities are usually far from exponential.Moreover, the impact of departures from exponential distributions on epidemic properties is usually greater in models incorporating small mixing groups, such as households.Although it is possible in principle to use the method of stages to extend the deterministic model in Sect.2.4 to include Erlang distributed infectious and latent periods, and to allow for varied local and global infection rates between stages of infection, in practice, the number of ODEs soon becomes infeasible.However, it is straightforward to calculate h D for such models, and indeed for models in which individuals have infectivity profiles (for example, Goldstein et al. (2009)), since such calculation only requires final outcome properties of an epidemic.We have found that h D > h D in most of our numerical studies with exponentially distributed infectious and latent periods, and that the difference is typically small unless the mean household size is large.We expect a similar conclusion to hold for models with other, more realistic, choices of infectious and latent period distributions.The different rows correspond to different choices of initial condition in terms of which household sizes contain initial infectives.In the first row, all initial infectives reside in size 1 households, and in the second row all initial infectives reside in size 2 households.In the third row, a fraction 2 households contain one initial infective and a fraction 2 households contain two initial infectives.As ↓ 0, h D ( ) approaches a value consistent with Table 2. Smaller values of give the same value as h D (10 −5 ) to this level of accuracy In Table 2 we present an example showing empirical evidence in support of this conjecture for the households SIR model when T E ∼ Exp(γ ).We also see in the bottom row of Table 2 evidence that mvar(H (m) D ) converges to a constant as m → ∞, which is consistent with a central limit theorem for H (m)  D .Recall from Sect.3.2.1 that we calculate h L D by considering a sequence of deterministic epidemics in which the initial fraction of the population that is infected tends to 0. Further, we conjecture that h L D ( (k) ) → h L D as k → ∞ for any sequence ( (k) ) satisfying (s,e,i,r )∈H (nmax) s,e,i,r ↓ 0 and (s,e,i,r )∈H 0rec

(k)
s,e,i,r ↑ 1 as k → ∞.We give an example for an SIR model in Table 3.
Tables 2 and 3 suggest that lim m→∞ H (m) D can be computed as the solution of the appropriate corresponding deterministic equations, along with the stopping condition R V (t) = 1. Figure 12 gives an example showing that, as m increases, the trajectories {S (m) (t)} and {R (m) V (t)} become smoother, with a random time translation capturing the time it takes for the stochastic epidemic to take off.As a result, the trajectories of R (m) V (t) and 1 − S (m) (t) against t can, for large m, be thought of as random time shifts of the corresponding deterministic trajectories.The random time shift is absent, however, when considering R     We compare R DI (π ) with the corresponding reproduction number R U (π ), when this fraction z of the population is vaccinated uniformly at random.Using (4.1) (cf.(4.2)), Then h n (π ) > 0 for π ∈ (0, π * n ( p)) and h n (π ) < 0 for π ∈ (π * n ( p), 1).Hence, h n (π ) = 1 − p has a unique solution, πn ( p) say, in (0, 1).Further, h n (π ) < 1 − p for π ∈ (0, πn ( p)), implying R U (π ) > R DI (π ), and h n (π ) > 1 − p for π ∈ ( πn ( p), 1), implying R U (π ) < R DI (π ).It follows that there exists λ
We are now in a position to prove Theorem 4.8.For ease exposition in the following proof, we denote h C (λ L ) by c(λ L ).Similarly, we denote h D (λ L ) by z(λ L ).Recall that n max < ∞, so all sums in the proof contain only finitely many terms and hence are easily differentiated.
Proof of Theorem 4.8 Suppose a fraction c(λ L ) are vaccinated prior to an epidemic, such that the epidemic becomes critical.Note c( 0 ).Now considering disease-induced herd immunity, assume that a fraction z(λ L ) are infected in a first epidemic such that the second epidemic is critical.Let π(λ L ) be the proportion who avoid global infection in the first epidemic.Then z( 0 from which Theorem 4.8 follows immediately by Taylor's theorem.We begin by considering the derivatives of c(λ L ) at λ L = 0.The post-vaccination threshold parameter RU (c) defined at (3.1) satisfies RU (c(λ L )) = 1, so by substituting i = n − v in the inner sum in (3.1) we have that Then n i=1 q n,i (c) = 1 − c, n i=1 q n,i (c) = −1 and n i=1 q n,i (c) = 0, by exchanging the order of derivative and summation.We also have Turning now to z(λ L ), let Pn,i (λ L , π(λ L )) = P( Sn = i) be the probability i members of a household avoid infection in the epidemic Ẽn (λ L , π(λ L )).Suppose all individuals infected (no longer susceptible) in the first epidemic are immune to infection in the second epidemic.Suppose the second epidemic is critical, so that R DI (π(λ L )) = 1.Then considering the remaining susceptibles yields We begin by making the notation in Theorem 4.11 explicit for the case n = 3.For i ∈ {0, 1, 2, 3}, P D i is the probability of i members being infected in a household of size 3 during an epidemic in the households model in which a proportion z are infected in the first epidemic.Similarly, P U i is the probability of a household containing i vaccinated individuals, assuming vaccination uniformly at random with probability z.Note that P D i , P U i and z are considered as functions of π ∈ (0, 1).We begin with a preliminary lemma.Write q 1 = φ(λ L ) and q 2 = φ(2λ L ) and observe that 0 < q 2 < q 1 < 1 for all λ L > 0. Further, P U 0 = (1 − z) 3 , P U 1 = 3z(1 − z) 2 and P U 2 = 3z 2 (1 − z).The system in (2.5), or direct calculation, gives P D 0 = π 3 , P D 1 = 3π 2 (1 − π)q 2 and P D 2 = 3π(1 − π)q 1 (2π(q 1 − q 2 ) + (1 − π)q 1 ).
We now prove Theorem 4.11 in the case n = 3.
Proof We show that RDI (z) > RU (z), from which the desired result follows.

C Early exponential growth rate r of SEIR households model
The early exponential growth rate (Malthusian parameter) r is defined as follows.Consider the branching process of infected households described in Sect.2.2.Then the Malthusian parameter, r , for the branching process satisfies where β(t) is the mean rate of global contacts emanating from a typical singlehousehold epidemic, t time units after the household is infected.For many choices of distributions for T E and T I , the left-hand-side of (C.1) is not tractable; in Ball et al. (2016) Section 2.8, ways to approximate it are considered.We restrict attention to the case T E ∼ Exp(δ) and T I ∼ Exp(γ ), in which case the dynamics are Markovian and the left-hand side of (C.1) can be computed as outlined below, cf.Pellis et al. (2011), Section 4.2, which considers an SIR model.
In the branching process in Sect.2.2, individuals correspond to infected households, and an gives birth whenever a global contact arises from the corresponding single-household epidemic.We extend Ball and Shaw (2015), Section 4.1, to include an exposed state.For n = 1, 2, . . ., let Ẽ(n) H denote a typical single-household epidemic in a household of size n, initiated by one member becoming exposed at time t = 0.For t ≥ 0, let S (n)  H (t), E (n) H (t) and I (n) H (t) denote respectively the number of susceptibles, exposed and infected individuals in Ẽ(n) H at time t.Let F (n) = {(s, e, i) ∈ Z 3 + : s + e + i ≤ n} be the set of possible household states for a household of size n.For (s, e, i) ∈ F (n) and t ≥ 0, let p   Let S n = F (n) = n 6 n 2 + 6n + 11 .For a given household of size n, we have, in an obvious notation, β (n) (t) = λ G (s,e,i)∈F (n)   i p s,e,i (t), which after conditioning on the size of a typical household in the approximating branching process yields αn (s,e,i)∈F (n)   i p s,e,i (t) dt.

Fig. 2
Fig. 2 Heat maps of the percentage error 100 h D − h D /h D as a function of ( p L , R 0 ) for common household sizes n = 2, 3, 4 and 5

Fig. 3
Fig.3Values of ĥ D with global restrictions scaled by a factor κ for the duration of the first epidemic, using the UK household size distribution (solid lines) and Morocco's household size distribution (dashed lines), taking λ L = 0.5 and considering several values of R 0 ) for π numerically and thereby determine λ * G (n, p), and the corresponding value R * 0 (n, p) of R 0 (which is independent of E[T I ]), such that the change in the ordering of h D and h C occurs. Figure 4 shows R * 0 (n, p) as a function of p for various n, with h D < h C above the plotted line and h D > h C below it.If p ≤ (n − 2)/2(n − 1) then no change of sign occurs and h D > h C for all values of R 0 > 1.

Fig. 8
Fig. 8 Critical values of (μ H , σ 2 H ) from Corollaries 4.5 and 4.9, together with values of these quantities for several countries' household size distributions

Fig. 9
Fig. 9 Herd immunity levels maintaining a fixed growth rate r implied by a given value of R BBT 0 , for the UK household size distribution

Fig. 10
Fig. 10 Herd immunity levels maintaining a fixed growth rate r implied by a given value of R BBT 0 , for Morocco's household size distribution

Fig. 11
Fig. 11 Comparison of h D , h D , h D − h D and h C respectively by country for R BBT 0 = 3 (with r held fixed) comparing Iraq, Chad, Morocco, UK, Japan and Finland

V
(t)  against 1 − S (m) (t), since the time shifts of R

Fig. 12
Fig. 12 Realised trajectories of R V (t) and 1 − S (m) (t) for m = 100 (left column) and m = 10000 (right column), with (λ G , λ L , γ ) = (2, 0.25, 1) and an equal number of households of size 1 and size 2 with a single initial infective residing in a household of size 2. The top row plots the trajectories of 1 − S (m) (t) against t.The second row plots the trajectories of R

V
(t) against t.The final row plots the trajectories of R

V
(t) against 1 − S (m) (t), with the epidemic progressing from left to right and with the corresponding deterministic trajectory superimposed 1 − S (m) (t) are the same for a given realisation of the process.In light of the definition above of H (m) D , we thus find, for large m, that H (m) D lim →0 h D ( ).Owing to the random time shift, however, it is not the case that T (m) * t * even for large m.B Proofs B.1 Proof of Theorem 4.7 Proof Let π ∈ (0, 1).Using (4.5), R DI (π ) = λ G E[T I ][(1 − p)π + npπ n ].Note that the proportion z infected in the first epidemic satisfies z = 1 − (1 − p)π − pπ n .

Table 1
The maximum percentage error for hD approximating h D when ( p L , R 0 ) ∈ [0, 1] × (1, 25],together with the parameter values for which the maximum is attained The following corollary, which involves only the mean and variance of H , is an immediate consequence of Theorem 4.4 in the case l * = 2. Corollary 4.5 If var( H ) < E[ H ]E[ H − 1] then h D > h C forhighly locally infectious epidemics which are just above threshold.If var( H ) > E[ H ]E[ H − 1] then h D < h C for highly locally infectious epidemics which are just above threshold.

Table 2
Empirical evidence the convergence of H D to a limiting value as m → ∞, using 10000 simulations of major outbreaks with (λ G , λ L , γ ) = (2, 0.25, 1) and an equal number of households of size 1 and size 2 The mean appears to converge toward a fixed value, with the standard deviation appearing to decrease toward 0 as m → ∞

Table 3
Empirical evidence of convergence of h D ( ) to h D , using (λ G , λ L , γ ) = (2, 0.25, 1) and an equal number of households of size 1 and size 2, with an initial fraction of infected households, where = (s,e,i,r )∈H(nmax)