An epidemic model with short-lived mixing groups

Almost all epidemic models make the assumption that infection is driven by the interaction between pairs of individuals, one of whom is infectious and the other of whom is susceptible. However, in society individuals mix in groups of varying sizes, at varying times, allowing one or more infectives to be in close contact with one or more susceptible individuals at a given point in time. In this paper we study the effect of mixing groups beyond pairs on the transmission of an infectious disease in an SIR (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}→ 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, both through a branching process approximation for the initial stages of an epidemic with few initial infectives and a functional central limit theorem for the trajectories of the numbers of infectives and susceptibles over time for epidemics with many initial infectives. We also derive central limit theorems for the final size of (i) an epidemic with many initial infectives and (ii) a major outbreak triggered by few initial infectives. We show that, for a given basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0, the distribution of the size of mixing groups has a significant impact on the probability and final size of a major epidemic outbreak. Moreover, the standard pair-based homogeneously mixing epidemic model is shown to represent the worst case scenario, with both the highest probability and the largest final size of a major epidemic. Supplementary Information The online version contains supplementary material available at 10.1007/s00285-022-01822-3.


Introduction
An assumption of the general SIR (susceptible → infective → recovered) epidemic model is that infection occurs through the interactions of pairs of individuals, one of whom is infectious and the other of whom is susceptible. The rate at which infections occur is governed by the rate at which pairs of individuals interact and the probability of transmission during an interaction. These two quantities are often subsumed into a single rate of infectious contacts made by an infective, with any infectious contact with a susceptible resulting in infection. This assumption of the epidemic being driven by pairs of interactions is, at least implicit, in most infectious disease models. The introduction of population heterogeneities in the form of variability in infectivity and susceptibility (eg. Neal 2007), metapopulation (eg. Ball and Clancy 1993), household (eg. Ball et al. 1997) or network (eg. Andersson 1998) models do not depart from the key assumption that infection is between pairs of individuals. One exception is the highly infectious household model, see eg. Becker (1995), where once an individual is infected in a household, their whole household becomes infected. However, this model is often framed in terms of allowing the within-household infection rate λ L → ∞, and as such represents the limit of the standard household model. Another exception is the Greenwood chain-binomial model, see eg. Bailey (1975), Chapter 14, which is a discrete-time model in which the probability that a given susceptible is infected at a given time depends only on their being infection in the population and not on the number of infectives. However, that assumption is appropriate only for small populations such as households.
The role of superspreaders or superspreading events is often mentioned in connection with the emergence of a disease, see eg. Lau et al. (2017) for Ebola 2014-2015 and Lewis (2021) for . Superspreaders are typically taken to be individuals who are particularly infectious rather than having atypical contacts, see eg. Lloyd-Smith et al. (2005), although in a network setting the infectiousness of an individual is linked to the number of neighbours they have in the network. In many countries one of the most significant NPI (non-pharmaceutical interventions) to reduce the spread of Covid-19 was the introduction of limits on the size of gatherings outside the home, see Flaxman et al. (2020). For example, in the UK, mass gatherings such as sporting events were prohibited and groups were limited in size to 6 individuals, with similar measures employed elsewhere. The motivation behind such an NPI is that gatherings of individuals allow the transmission of a disease from a single infective or small group of infectives to a large group of susceptible individuals. Mass gatherings, such as a football match, typically take place over a short time period, say two hours, in comparison to the infectious period of an individual which will usually be several days. Moreover, individuals infected at a gathering taking place over a time period of two hours are unlikely to become infectious and start transmitting the disease during the same gathering. Therefore, we extend the standard SIR epidemic model to assume that individuals, rather than making contact with single individuals at the points of a Poisson process, are involved in mixing events at the points of a Poisson process with mixing events involving c ≥ 2 individuals. If there is at least one infective and at least one susceptible in a mixing event then there is the possibility of an infection taking place but multiple infections can take place within a single mixing event. It is necessary to model the rate at which mixing events occur and the distribution of the size of mixing events separately from the transmission of the disease within mixing events.
In this paper we consider both the early stages of the epidemic through a branching process approximation and the trajectory of infectives and susceptibles through time in the event of a major epidemic outbreak. By studying the epidemic process as n → ∞, we are able to compute key epidemic quantities such as the basic reproduction number, R 0 , and the probability of a major outbreak from the branching process approximation, and the limiting normal distribution for the final size of a major outbreak through a time transformed process. We show that for a given R 0 , both the probability and mean final size of a major outbreak depend on the distribution of the sizes of mixing events, with the standard SIR epidemic model with pairwise interactions having both the highest probability and the largest mean size of a major outbreak amongst all mixing event distributions. Moreover, for a given R 0 > 1 and a single initial infective, the probability of a major outbreak can vary between 0 and 1 − R −1 0 (its value for the standard SIR epidemic model) and the mean final size (fraction of the population infected) can vary between 1 − R −1 0 and τ * , where τ * solves 1 − τ * = exp(−R 0 τ * ) and is the mean final size of the standard SIR epidemic model. For example, if R 0 = 2, the probability of a major outbreak can range between 0 and 0.5 and the mean final size of a major outbreak can range between 0.5 and 0.7968.
The remainder of the paper is structured as follows. In Sect. 2, we define the SIR epidemic model with mixing events. The main results are summarised in Sect. 3, which include the probability of a major outbreak along with the mean and variance of the size of a major outbreak for a general mixing event distribution. These results are asymptotic as the population size n → ∞. In Sect. 4, we obtain important inequalities which allow us to provide more generally lower and upper bounds for the probability of extinction and the mean final size, along with orderings of mixing event distributions enabling us to show the above-mentioned comparisons with the standard SIR epidemic model. In Sect. 5, we consider logarithmic and geometric mixing event distributions which are amenable to analysis and allow us to obtain explicit expressions for the extinction probability, which further highlight the role of the mixing event distribution. In Sect. 6, we illustrate numerically a selection of the results and use simulations to show their relevance for finite population sizes, n. The proofs of the results presented in Sect. 3 are given in Sect. 7. In Sect. 8, we give some concluding comments and discuss possible directions for future work.

Model
We consider an SIR epidemic in a population of n individuals in which infection occurs via mixing events that occur at the points of a Poisson process having rate nλ. For i = 1, 2, . . . , the i th mixing events involves C  infective present has probability π c of making an infectious contact with any given susceptible present, with all such contacts occurring independently. Mixing events are assumed to be instantaneous. Any susceptible who is contacted by at least one infective at a mixing event becomes infected and remains so for a time that follows an Exp(γ ) distribution, i.e. an exponential distribution with rate γ and hence mean γ −1 . There is no latent period but newly infected individuals cannot infect susceptibles during the mixing event at which they were infected. All processes and random variables involved in the above model are mutually independent. The process starts with m n individuals being infective at time t = 0, with the remaining n − m n being susceptible, and ends when there is no infective remaining in the population. Denote this epidemic model by E (n) .
An example of a realisation of an epidemic outbreak is presented in Fig. 1 with n = 10, m n = 1, P(C (n) = 3) = 1 and π 3 = 1, so all mixing events are of size 3 and all susceptibles at a mixing event that contains at least one infective are necessarily infected at that event. There are 7 mixing events, labelled (a) to (g) in chronological order, during the time period shown. Mixing event (a) contains only susceptibles so no new infectives occur. Mixing event (b) involves the initial infective, individual 1, and two susceptibles, individuals 2 and 4, so the latter two become infected. Mixing event (c) involves two infectives, individuals 2 and 4, and one susceptible, individual 7, so individual 7 becomes infected. Mixing event (d) involves one infective, individual 2, one recovered, individual 4, and one susceptible, individual 5, so the latter becomes infected. Mixing event (e) involves one infective, individual 5, and two recovered, individuals 4 and 7, so no new infectives occur. The only remaining infective, individual 5, recovers between mixing events (e) and (f), so the epidemic terminates.
Note that if P(C (n) = 2) = 1, so all mixing groups necessarily consist of 2 individuals, then E (n) is the standard homogeneously mixing stochastic SIR epidemic model, with recovery rate γ and individual-to-individual infection rate 2λπ 2 n−1 . (If there are s susceptibles and i infectives in the population, the probability that a mixing group of size 2 contains one infective and one susceptible is si/ n 2 = 2si n(n−1) , so the rate at which new infections occur is nλ × 2si n(n−1) × π 2 = 2λπ 2 n−1 si.) Note also that the groups formed at mixing events are very different from the small mixing groups in for example the households model. In the latter, these groups are permanent, while in the present model their duration is instantaneous.

Introduction
We now present the main results of the paper, which are concerned with the behaviour of the epidemic model introduced in Sect. 2 as the population size n → ∞. We assume that C (n) D −→ C as n → ∞, where D −→ denotes convergence in distribution and C has probability mass function p C (c) = P(C = c) (c = 2, 3, . . . ). In Sect. 3.2 we consider epidemics with many initial infectives, more precisely the asymptotic regime n −1 m n → as n → ∞, where > 0. Thus in the limit as n → ∞ a strictly positive fraction of the population is initially infected. We give a law of large numbers (LLN), which puts a deterministic approximation of the model on a rigorous footing (Theorem 3.1) and an associated functional central limit theorem (CLT), which describes fluctuations about the limiting deterministic model for large n (Theorem 3.2). We also give a CLT for the final size of an epidemic (Theorem 3.3).
In Sect. 3.3 we consider epidemics with few initial infectives, more precisely the asymptotic regime in which m n = m for all sufficiently large n. Thus, for such n, each epidemic has m initial infectives. Under these conditions, for suitably large n, the early stages of the epidemic E (n) can be approximated by a branching process, which we use to determine R 0 and the probability that, in the limit n → ∞, an epidemic becomes established and leads to a major outbreak which infects a strictly positive fraction of the population. We also present a CLT for the size of a major outbreak (Theorem 3.5). The conditions for some of the theorems are quite involved, so in Sect. 3.4 we give simple, easily checkable sufficient conditions which cover most, if not all, practical applications.

Epidemics with many initial infectives
Throughout this section we assume that n −1 m n → as n → ∞, where > 0.

Temporal behaviour
For t ≥ 0, let S (n) (t) and I (n) (t) be respectively the numbers of susceptibles and infectives in E (n) at time t. Throughout the paper, for a vector x ∈ R d for some specified d, |x| denotes its Euclidean norm. Further, p −→ denote convergence in probability.
where (x(t), y(t)) is given by the solution of the following system of ODEs (ordinary differential equations) with initial condition (x(0), y(0)) = (1 − , ): The ODEs (3.1) yield a deterministic approximation to the epidemic E (n) . The next result is concerned with fluctuations of the sample paths of E (n) about that deterministic limit. Let where I denotes the 2×2 identity matrix. For n = 2, 3, . . . , let p (n) where (x(t), y(t)) is as in Theorem 3.1. Then where ⇒ denotes weak convergence in the space of right-continuous functions f : [0, ∞) → R 2 having limits from the left (i.e. càdlàg functions), endowed with the Skorohod metric, and {V (t) : t ≥ 0} is a zero-mean Gaussian process with V (0) = (− 0 , 0 ) and covariance function given by where denotes transpose.
Note that it follows immediately from (3.12) that Differentiating (3.13) and using (3.9) yields that (t) satisfies the matrix ODE with initial condition (0) = 0, the 2 × 2 matrix of zeros. Thus (t) can be computed by numerically solving the ODEs (3.1) and (3.14) simultaneously.

Final outcome
Let T (n) = n−S (n) (∞) be the total number of individuals infected during the epidemic (including the initial infectives), i.e. its final size. As detailed in Sect. 7.4, we prove a CLT for T (n) by considering a random time-scale transformation of {(S (n) (t), I (n) (t)) : t ≥ 0} in which the clock is slowed down by a factor n −1 I (n) (t). This leads to the time-transformed deterministic model given by (3.16) Adding the differential equations (3.15) and solving yields (3.18) The time transformation does not change the final outcome of the epidemic. Let τ = inf{t > 0 :ỹ(t) = 0}, so using (3.17),x(τ ) = 1 − γτ . Note that x(∞) = x(τ ), so the fraction of the population infected (including the initial infectives) in the deterministic epidemic given by (3.1), with (x(0), y(0)) = (1 − , ), is γτ . The CLT for T (n) is proved by deriving an analogous functional CLT to Theorem 3.2 for the time-changed process and then considering an associated crossing problem. Before stating the theorem some more notation is required. Let where R 0 , the basic reproduction number (see Sect. 3.3), is given by (3.20) Let (σ 2 S, (t),σ S I , (t),σ 2 I , (t)) be the solution of the following system of ODEs with initial condition (σ 2 S (0),σ S I (0),σ 2 I (0)) = (0, 0, 0): (3.21) where (x(t),ỹ(t)) are given by (3.17) and (3.18), and For future reference, defineτ 0 and σ 2 T (0) by setting = 0 in the above. For i, j = 0, 1, . . . and y 0 > 0, let The CLT assumes that the following conditions hold: and there exists y 0 > 0 is such that where N(0, σ 2 T ( )) denotes a normal distribution with mean 0 and variance σ 2 T ( ).

Epidemics with few initial infectives
Throughout this section, we assume that m n = m for all sufficiently large n. Consider a typical infective, i * say, in the early phase of the epidemic E (n) . The probability that a mixing event of size c involves i * is c n . Mixing events occur at rate nλ and have sizes that are i.i.d. realisations of C (n) , so mixing events that involve i * occur at rate λμ , and, for c = 2, 3, . . . , n, the probability that a mixing event is of size c given that it includes i * is cp C . Moreover, since mixing events are formed by choosing individuals uniformly at random from the population, with probability close to 1, every mixing event involving i * consists otherwise entirely of susceptibles. By considering the limits of the above quantities as n → ∞, the early phase of the epidemic E (n) can be approximated by the following branching process, denoted by B, in which individuals and birth events correspond respectively to infectives and mixing events in the epidemic process E (n) .
Individuals in the branching process B behave independently. A given individual, i * say, has lifetime L ∼ Exp(γ ), during which they have birth events at the points of a Poisson process having rate λμ C . The numbers of offspring individual i * produces at their successive birth events,Z 1 ,Z 2 , . . . , are i.i.d. copies of a random variableZ having distribution defined as follows. The size of a typical mixing event involving i * is distributed asC, the size-biased version of C, having probability mass function If the mixing event is of size c then there are c − 1 susceptibles present, each of whom is infected independently with probability π c by i * . Thus,Z is distributed as a mixture of Bin(c − 1, π c ) (c = 2, 3, . . . ) random variables with respective mixing probabilities pC (c) (c = 2, 3, . . . ), Note that an individual may have no offspring at a birth event.
Let R denote the total number of offspring a typical individual has in B. Then where G is the number of birth events that an individual has in their lifetime and R = 0 if G = 0. Further, G is independent ofZ 1 ,Z 2 , . . . and has the geometric distribution (3.33) Hence, where R 0 is given by (3.20). Hence, as stated in Sect. 3.2.2, R 0 is the basic reproduction number of the epidemic E (n) . Let z denote the extinction probability of B given that initially there is one individual. Then from standard branching process theory, z is given by the smallest solution in .) It follows from (3.32) and (3.33) that (3.36) The Malthusian parameter (exponential growth rate) r of the branching process B is easily obtained since the mean rate an individual produces offspring t time units after their birth is P( Note that if R 0 and γ are held fixed, then r is the same for all corresponding choices of the distribution of C and π k (k = 2, 3, . . . ). In particular, under these conditions, the early exponential growth of an epidemic that takes off is the same as that of a standard homogeneously mixing epidemic.
The formulae for R 0 and f R (s) simplify when the infection probability π c is independent of c, say π c = π for all c. In particular, (3.34) then yields The approximating branching process B can be put on a rigorous footing by theorems concerning convergence of E (n) to B as n → ∞. However, their proofs are long and are presented separately in Ball and Neal (under review). The usual approach to proving such theorems (see e.g. Ball and Donnelly 1995) is to construct sample paths of E (n) for each n from those of the limiting branching process B. That approach is not easily implemented in the present setting as mixing groups induce dependencies between infectives. The following theorem, that is used in the proof of Theorem 3.5 below, is proved in Ball and Neal (under review). and conditions (3.10) and (3.28) hold. Then

Theorem 3.4 (a) Suppose that m n = m for all sufficiently large n, C
In view of Theorem 3.4(a), we define a major outbreak to be one whose final size is at least log n. Theorem 3.4(b) implies that a major outbreak infects at least a fraction δ of the population with probability tending to one as n → ∞. However, δ depends on the parameters of the epidemic E (n) and can be arbitrarily close to 0. The final result in this section gives a CLT for the size of a major outbreak. Recall the notationτ 0 and σ 2 T (0) introduced just after (3.24).

Sufficient conditions for theorems
Simpler, easily checkable sufficient conditions for the above theorems are proved in Appendix A in the Supplementary Information. If there is a maximum mixing group size, c * say, so P(C (n) > c * ) = 0 for all n, then the conditions concerning C (n) and C of all of the theorems are satisfied if lim n→∞ √ n[ p . Two natural choices for the distribution of C (n) when there is no maximum mixing group size are C (n) D = min(C, n) and C (n) D = (C|C ≤ n), where D = denotes equal in distribution. In both of these cases the conditions concerning C (n) and C are satisfied if f C (s 1 ) < ∞ for some s 1 > 1. This condition requires that C has finite moments of all orders but it holds for all choices of infection probabilities π c (c = 2, 3, . . . ). If these probabilities satisfy π c ≤ ζ c for all sufficiently large c, for some ζ ∈ (0, ∞), then E[C 4 ] < ∞ is sufficient for the conditions concerning C (n) and C.

Constant mixing event infection probability c
In this subsection we assume that π c = π for all c = 2, 3, . . . . Recall that z denotes the extinction probability of the branching process B. Further, for ∈ (0, 1], for a sequence of epidemics (E (n) ) satisfying lim n→∞ n −1 m n = ,τ is the fraction of the population that is infected by the epidemic E (n) in the limit as n → ∞. If instead, m n = m for all sufficiently large n, thenτ 0 is the fraction of the population that is infected by a major outbreak in the limit as n → ∞.
The following proposition facilitates comparison of the above properties between epidemics when R 0 is held fixed and also calculation/computation of some of these properties in the special cases considered in Sect. 5. Suppose that E[C 2 ] < ∞ and let C be a random variable having probability mass function Assume without loss of generality that the time scale is chosen so that γ = 1, The same argument as in the proof of part (a) now shows that λg(ỹ) = 1 y U (ỹ), as required.
Let D denote the space of random variables taking values in {2, 3, . . . } and, for c = 3, 4, . . . , let D c be the subspace of D consisting of random variables that take let Const(c) denote the distribution consisting of a unit mass at c.
Suppose that C ∼ Const(2). Then U (y) = R 0 y (0 ≤ y ≤ 1) and it follows using Of course these results coincide with those of the standard homogeneously mixing SIR epidemic model. Note that Proposition 4.1 implies that z andτ ( ∈ [0, 1)) are determined by Proof (a) Write U (y) as U (y, C, π) to show explicitly its dependence on (C, π).
Theorem 4.1(a) shows that if R 0 > 1 and the event size distribution C are held fixed, then the epidemic increases as the infection probability π decreases, in the sense that both the probability of and the fraction of the population infected by a major outbreak increase. Moreover, as π ↓ 0, these properties tend to the corresponding properties of the standard homogeneously mixing SIR model. (If C ∼ Const (2), so all mixing events are of size 2, and R 0 is held fixed, then the distribution of the final outcome of the epidemic is independent of π ∈ (0, 1].) Theorem 4.1(b) shows that if R 0 > 1 and π ∈ (0, 1] are held fixed, then the epidemic tends to decrease (in the above sense) as the mixing event size C increases. Moreover, the standard homogeneously mixing SIR model necessarily provides an upper bound for the model with mixing events. Furthermore, if there is a maximum mixing event size, c max say, then the model in which mixing events have size c max yields a lower bound.
is the size (including the initial infectives) of the standard deterministic homogeneously mixing SIR epidemic having basic reproduction number R 0 , when the initial fractions susceptible and infected are R −1 0 and , respectively. Let P denote the set of all possible π.

Theorem 4.2 (a) For any
Then, for any (C, π ) ∈ D × P with R 0 > 1, and Proof Throughout the proof we assume without loss of generality that γ = 1.
Turning to (4.7), note that fZ (s) > 0 for s ≥ 0, since R 0 > 1 precludes P(Z = 0) = 1. Hence, recalling λ = a(C, π )R 0 /μ C , it follows from (3.35) that f R (s) > 1/[1 + a(C, π )R 0 ] for s ∈ [0, 1]. The lower bound for z in (4.7) follows, since z is the smallest solution in [0, 1] of f R (s) = s. The upper bound is immediate as R 0 > 1. Theorem 4.2(a) shows that if R 0 > 1 is held fixed, the standard homogeneously mixing epidemic (C ∼ Const (2)) always provides an upper bound for the model with mixing events. Theorem 4.2(b) provides useful bounds when mixing group sizes can be large. Suppose that C ∼ Const(c) and π c = π ∈ (0, 1]. Then a(C, π ) = 1 (c−1)π . It follows that when c is large, the probability of a major outbreak is small (0 in the limit c → ∞) and if one occurs its size is slightly larger than 1 − 1 R 0 (its limit as c → ∞). Thus for large c, in the event of a major outbreak, the epidemic effectively stops spreading once it becomes critical. More generally, these conclusions hold if a(C, π ) is small.
The proof of Theorem 4.2(b) also has implications for the duration of an epidemic. The duration of the standard homogeneously mixing epidemic is studied in Barbour (1975). For epidemics with few initial infectives, a major outbreak can be split into three phases: an initial phase, the main body of the epidemic and a final phase, during which the cumulative number infected increases to 1 n, (τ − 2 )n and τ n, respectively, where 1 and 2 are both small. For large n, the durations of the initial and final phases are both of exact order (log n), and the duration of the middle phase is of exact order (1), more precisely the time for x(t) to decrease from 1 − 1 to τ + 2 in the limiting deterministic model. A similar decomposition holds for the model with mixing events. The durations of the initial and final phases of a major outbreak, during which the epidemic can be approximated by a supercritical and a subcritical branching process, respectively, are again of exact order (log n). The proof of Theorem 4.2(b) implies that y(t) ≤ a(C, π )R 0 throughout the middle phase, so the duration of the middle phase is very long if a(C, π ) is small, and tends to infinity as c → ∞ in the above case when all mixing events have size c.

Logarithmic event distribution
Suppose that C follows a logarithmic distribution with parameter α ∈ (0, 1), i.e. (5.1) Hence, using (4.2), and application of Proposition 4.1(a)(ii) yields that, when R 0 > 1, Now fĈ −2 (s) is increasing in α for any s ∈ [0, 1] and strictly increasing for s ∈ [0, 1). It follows using Theorem 4.1(b)(i) that if R 0 > 1 and π are held fixed then both the probability and final size of a major outbreak are increasing in α. Moreover, these epidemic properties tend to those of the standard SIR model as α → 1. Note also that, in the notation of Theorem 4.2(b), a(C, π ) = α π . Thus, by that theorem, the fraction of the population infected by a major outbreak tends to 1 − R −1 0 as α ↓ 0.

Geometric event distribution
Suppose that C follows a Geometric distribution with parameter α ∈ (0, 1], conditioned to be strictly greater than 1, so . Hence, using (4.2), Let w be the largest solution in [0, 1] of (1 − y)U (y) = y. A little algebra shows that, if R 0 > 1, then w satisfies a quadratic equation and application of Proposition 4.1 yields The probability-generating function fĈ −2 (s) is increasing in α for any s ∈ [0, 1], so the behaviour of the epidemic properties with α, noted in Sect. 5.1 for logarithmic event distributions, hold also for geometric event distributions. Note that, in an obvious notation, for fixed α ∈ (0, 1), f s) (0 ≤ s < 1). Hence, by Theorem 4.1(b)(i), if R 0 > 1 and α are held fixed, then both the probability and final size of a major outbreak are greater when the sizes of mixing events follow a logarithmic distribution than when they follow a geometric distribution.

Numerical illustrations
In this section we present numerical results which demonstrate the usefulness of our CLTs for finite population size n and illustrate the mixing group size distribution, mean mixing group size and π affect model properties. The recovery rate γ = 1 throughout this section.
In Fig. 2, we consider epidemics in a population of size n = 100, 000 with examples of a logarithmic and a geometric mixing event distribution with R 0 = 2. For the two mixing event distributions, we plot the trajectories of infectives for the time interval [0, 20], given by the solution of the ODE (3.1) scaled by the population size n, together with 95% equal-tailed probability intervals obtained using the functional CLT, Theorem 3.2 with = 0.001 (100 initial infectives). Superimposed on each mixing event distribution plot are 100 simulations of the epidemic with 100 initial infectives on the interval [0, 20]. There is very good agreement between the simulations and the asymptotic limits with excellent coverage of the probability intervals by the simulated trajectories. We note that for the geometric mixing event distribution two of the simulations experience a delay in becoming established and consequently their trajectories have peak later. This is a consequence of variability in the early stages of an epidemic before the number of infectives has increased to a level at which the deterministic approximation fully takes hold.
In Table 1, we show how as n → ∞ the mean proportion (n −1 E[T (n) ]) and scaled variance (n −1 var(T (n) )) of the final size of a major outbreak approach their asymptotic limits obtained from Theorem 3.5. For fixed size mixing groups of size C = 3, 10,000 major outbreaks, each initiated by a single infective, were simulated for each of the following combinations of R 0 (= 1.2, 1.5, 2.0, 3.0), π (= 0.2, 0.5, 1.0) and n (= 10 4 , 10 5 , 10 6 ). For each set of parameter values, a cut-off for what constitutes a major outbreak for finite n is required. Figure 3 shows histograms of the proportion of the population ever infected in 10,000 simulated epidemics for each of n = 10, 000 and n = 100, 000, when R 0 = 1.2, π = 1 and C = 3, with each simulated epidemic being initiated by 10 infectives. (We initiated the epidemics with 10 rather than one infective, since with one initial infective very small minor outbreaks dominate the histograms.) There is a clear distinction between minor and major outbreaks when n = 100, 000 but when n = 10, 000 there is no clear distinction and the choice of cut-off is more arbitrary. The vertical red lines show the cut-off used to produce the corresponding entries in Table 1. The cut-offs for other choices of parameter values and n were determined by similar examination of histograms. More generally, the lack of distinction between minor and major outbreaks becomes more pronounced Table 1 Mean proportion infected and scaled variance for epidemics with mixing groups C   Proportion of the population infected in a major outbreak, n −1 T (n) . Histogram for 1,000 simulations in population of size n = 10 5 with R 0 = 2, π = 0.50 and a logarithmic mixing event distribution with α = 0.05 (μ C = 8.82). Superimposed is the probability density function of N (0.5972, 0.01184 2 ), the approximating normal distribution. See text for further details as the size of mixing groups increases for R 0 and π close to 1. Note from Table 1 that there is excellent agreement between the asymptotic and finite n values for the mean proportion of the population infected by a major outbreak. There is also good agreement for the variances for all choices of n unless R 0 is close to one, in which case n needs to be larger for the asymptotic results to provide a good approximation.
In Fig. 4, we illustrate the CLT for the final size using 1,000 simulations of major outbreaks, with a cut-off of = 1, 000, initiated by a single infective in a population of size n = 100, 000, with R 0 = 2, π = 0.50 and a logarithmic mixing event distribution with α = 0.05 (μ C = 8.82). The asymptotic values of n −1 E[T (n) ] and n −1 var(T (n) ) areτ 0 = 0.5972 and σ 2 T (0) = 14.0274, respectively, giving an approximate normal distribution for the proportion infected as N (0.5972, 0.01184 2 ). The superimposed normal approximation closely matches the histogram of the proportions infected in a major outbreak.
In Figs. 5 and 6 we show the trajectories of the mean proportion and scaled variance of susceptibles and infectives over the time interval [0, 20] for logarithmic mixing  Fig. 5, we vary π between 0.01 and 1 keeping α = 0.35 fixed, so μ C = 3.02. In Fig. 6, we fix π = 0.5 and vary α between 0.1 and 1, corresponding to mean mixing event sizes ranging from 5.78 down to 2 (C D −→ Const(2) as α ↑ 1). We observe that as the final size decreases the epidemic has a smaller peak and heavier tail corresponding to a longer duration of the epidemic. This is accompanied by greater variability in the number of susceptibles and infectives over time. Further, if we vary , the initial proportion infected, say to = 0.001, the effect on the mean trajectories is small, approximately corresponding to a time-shift in the trajectory, but the variance trajectories change significantly, with an increase in magnitude by a factor of between 5 and 10 at the peak and also a larger distinction between the two modes in the variance.
In Fig. 7, we demonstrate how the final size (asymptotic mean proportion infected) of the epidemic,τ 0 , varies with mean mixing event size, μ C , mixing event distribution (logarithmic, geometric and fixed size) and constant π for a fixed R 0 = 1.5. We observe that for all mixing event distributions and choice of π , the final size agrees for μ C = 2 (all events of size 2) andτ 0 decreases as μ C increases. The decrease inτ 0 as μ C increases is, as we would expect, more marked the larger π is. We also note that the logarithmic distribution has the fastest drop-off, followed by the geometric distribution, with the fixed group size experiencing the slowest reduction as the mean event size grows. We observe that for the logarithmic distribution the final sizeτ 0 < 0.35 for μ C = 20, approaching the lower bound of 1 − 1/R 0 = 0.3333.

Introduction
We prove these theorems by using the theory of density dependent population processes (see eg. Ethier andKurtz (1986), Chapter 11, andPollett (1990)). In Sect. 7.2, we explain how our model fits into that framework and state the general LLN and functional CLT that we use (Theorems 7.1 and 7.2, respectively). In Sect. 7.3, we derive the Fig. 7 Proportion of individuals infected,τ 0 , against the mean mixing event size, μ C , for R 0 = 1.5. Solid linelogarithmic distribution, dashed line -geometric distribution and points -fixed size distribution. π = 0.01 (black), 0.10 (red), 0.25 (green), 0.50 (blue) and 1.00 (purple) (color figure online) key limiting drift and variance/covariance functions for application of those theorems and prove that the conditions of these theorems are satisfied in their application to the temporal behaviour of epidemics with many initial infectives (Theorems 3.1 and 3.2). In Sect. 7.4, we consider the CLTs for the final size of epidemics (Theorems 3.3 and 3.5).

Theorem 7.1 Suppose that F is Lipschitz continuous on E,
Then, for any t 0 > 0, For (x, y) ∈ E define the infinitesimal variance/covariance matrices (n) l lβ (n) l (x, y) and G(x, y) = l∈ l lβ l (x, y), (7.7) and let The following theorem follows from Pollett (1990), Theorem 3.2. The covariance function for the limiting Gaussian process is taken from Ethier and Kurtz (1986), Chapter 11, equations (2.19) and (2.21).

Proofs of Theorems 3.1 and 3.2
We first determine the limiting drift and variance/covariance functions F and G. Recall that in the epidemic in a population of size n, the individuals involved in a given mixing event of size c are chosen by sampling c individuals uniformly at random from the population without replacement. In the limit as n → ∞ this converges to a corresponding sampling with replacement. For (x, y) ∈ E and c = 2, 3, . . . , let μ c (x, y) = E[Z ] and μ c,2 (x, y) = E[Z 2 ], where Z is the number of new infectives created in a mixing event of size c in which each of the c individuals is independently susceptible, infective or recovered with probabilities x, y and 1 − x − y, respectively.
We show via a sequence of lemmas that the conditions of Theorems 7.1 and 7.2 are satisfied. For i, j = 0, 1, . . . , let
(c) Note that 0 ≤ h(x, y) ≤ E[C 2 ] for all (x, y) ∈ E, so G is bounded on E. Omitting the details, elementary calculus yields that for all (x, y) ∈ E, Since E[C 2 ] < ∞ and S(3, 1) < ∞, it follows that the partial derivatives of h exist and are bounded on E, so h is uniformly continuous on E, whence so is G.
We consider now the convergence of F (n) and G (n) to F and G as n → ∞. Let E n = n −1 E (n) and note from (7.1) that the transition intensities of {(S (n) (t), I (n) (t))} require only the values of β (n) l (x, y) for (x, y) ∈ E n . Thus in defining the functions β (n) l (l ∈ (n) ) we are free to interpolate continuously between the points of E n rather than using the functional form arising from the appropriate hypergeometric distribution. Recall μ c (x, y) and μ c,2 (x, y) defined just before Lemma 7.1. For c = 2, 3, . . . , n = c, c + 1, . . . and (x, y) ∈ E n consider a mixing event of size c in a population of size n that contains s = nx susceptibles and i = ny infectives. Let Z n be the number of new infectives created at this mixing event, μ (n) Lemma 7.3 (a) For c = 2, 3, . . . , n = c, c + 1, . . . and all (x, y) can be defined using continuous interpolation between the points of E n , so that for c = 2, 3, . . . , n = c, c + 1, . . . and all (x, y) ∈ E, c are distinct. If D n occurs the spread of infection is coupled so that Z n = Z . The precise construction of (Z n , Z ) when D n does not occur is not needed for our argument. For an event, D say, 1 D denotes its indicator function. Now Also, The first inequality in (7.16) follows. The second one is proved similarly using 0 ≤ Z 2 , Z 2 n ≤ c 2 . (b) Fix n. Since the two inequalities in part (a) hold for all (x, y) ∈ E n , and β l (l ∈ ), μ c and μ c,2 are continuous on E, β (n) l (l ∈ (n) ) may be interpolated continuously between the points of E n so that the two inequalities in (7.17) hold for all (x, y) ∈ E.
For n = 2, 3, . . . and (x, y) ∈ E, let Observe, using (7.2) and (7.7), that (7.24) where at the final step we have used μ (n) c (x, y) ≤ c for all (x, y) ∈ E in the first inequality and Lemma 7.3(b) in the second inequality. Now lim n→∞ cp (7.25) As n → ∞, the first term on the right-hand side of (7.25) tends to 0 by assumption and the second term tends to 0 as E[C 3 ] < ∞, so (7.23) follows.
We now use the above lemmas to show that the conditions of Theorems 7.1 and 7.2 are satisfied. Note from the proofs of those theorems that the conditions on F, G, F (n) and G (n) only need to be satisfied in an open neighbourhood of the limiting trajectory ((x(t), y(t)) : t ≥ 0); see eg. the first line of the proof of Ethier and Kurtz (1986), Theorem 11.2.1. For any t 0 > 0, the trajectory ((x(t), y(t)) : 0 ≤ t ≤ t 0 ) is contained in the interior E o of E, since > 0, so we may use E o in the definition of {(S (n) (t), I (n) (t))} (n = 1, 2, . . . ) as an asymptotically density dependent family.
To prove Theorem 3.1, note that F is Lipschitz continuous on E by Lemma 7.2(a), since S(2, 1) < E[C 2 ] < ∞, and that condition (7.4) is satisfied by assumption. Further, recalling (7.20), Lemma 7.4(a) shows that condition (7.5) is satisfied. Thus all the conditions of Theorem 7.1 are satisfied and Theorem 3.1 follows.

Proofs of Theorems 3.3 and 3.5
The final size T (n) of E (n) is given by T (n) = n − S (n) (τ (n) ), where τ (n) = inf{t > 0 : I (n) (t) = 0}. We study the asymptotic distribution of T (n) via the exit of the process {(S (n) (t), I (n) (t))} from E (n) + , the set of states in E (n) with i > 0. Under the conditions of Theorem 3.3, τ (n) p −→ ∞ as n → ∞, so to obtain a process which exits E (n) + in finite time in the limit as n → ∞, we consider {(S (n) (t),Ĩ (n) (t))}, a random time-scale transformation of {(S (n) (t), I (n) (t))} in which, at any time t ≥ 0, the clock is slowed down by a factor n −1 I (n) (t). Thus {(S (n) (t),Ĩ (n) (t))} is a continuous-time Markov chain with transition intensities (cf. (7.1)) satisfying, for (s, i) ∈ E (n) + and l ∈ (n) , The use of such random time-scale transformations in epidemic modelling goes back to Watson (1980).
( 7.27) It is possible that β l (x, y) < 0 for some (x, y) with y < 0. If so, replace l by l = −l and define β l (x, y) = −β l (x, y) for such (x, y). This change, which is assumed implicitly in the following, does not effect the resulting infinitesimal drift and variance/covariance functions. In (7.27), y −1 β l (x, y) is defined by continuity when y = 0. It is immediate that the limit exists for recovery jumps, while for infection jumps the limit exists as for a mixing event to yield new infectives at least one of its members must be an infective. For l ∈ , defineβ l byβ l (x, y) = y −1 β l (x, y) ((x, y) ∈Ẽ(y 0 )).
The following theorem is proved via an analogous sequence of lemmas to those used in the proof of Theorem 3.2. The proofs of the lemmas corresponding to Lemmas 7.2 and 7.4 are obvious generalisations, although the details are more involved. The proof of the lemma corresponding to Lemma 7.3 is less immediate. Details may be found in Appendix B in the Supplementary Information.

Concluding comments
We have presented a new class of SIR epidemic models in which disease is transmitted via mixing groups, and not just pairwise interactions, together with a branching process approximation for the early phase of an epidemic with few initial infectives and CLTs for both the temporal behaviour of epidemics with many initial infectives and the final outcome of a major outbreak. The standard homogeneously mixing SIR epidemic model is a special case of the mixing group model when mixing groups are necessarily of size 2. If R 0 and the recovery rate γ are held fixed then the initial exponential growth rate r of our model is the same as that of the standard SIR model but both the probability, 1 − z, and final size, τ , of a major outbreak are smaller. We have proved a number of comparison results for the mixing group model. In broad terms, if R 0 and γ are held fixed, then both 1 − z and τ decrease with mixing group size and the duration of an epidemic increases with mixing group size. In the extreme case of very large mixing groups, the final size τ of a major outbreak can be only fractionally larger than 1 − R −1 0 . (Note that τ is necessarily greater than 1 − R −1 0 since otherwise the epidemic would never become subcritical.) One limitation of the analysis and results presented is that they assume that the infectious period follows an exponential distribution. The approximating branching process B described in Sect. 3.3 can be extended in the obvious fashion to the case when the infectious period follows an arbitrary but specified distribution. The CLTs presented in this paper can be extended via the method of stages to incorporate an Erlang, and more generally a phase-type, infectious period distribution. An alternative approach is to use the LLN and functional CLT for age and density dependent population processes in Wang (1975Wang ( , 1977, though details are likely to be complicated. The phase-type approach is generally easier to implement as it requires just the methodology of density dependent population processes, which is well developed. However, if the infectious period distribution does not belong to the class of phase-type distributions then that approach is only approximate since one needs to approximate the actual infectious period distribution by one which is phase-type. Although any distribution can be approximated arbitrarily closely by a phase-type distribution, the size of the phase space can become infeasible for numerical purposes. It would be interesting to study other models for transmission within mixing events. One possibility is to make a Greenwood-type assumption: viz. if there are c individuals at a mixing event then, provided there is at least one infective at the event, each susceptible present becomes infected independently with probability π c . The approximating branching process B is unchanged, as it assumes that all mixing events contain at most one infective. However, the CLTs do change. Omitting the details, the formulae for g c (y) and h c (x, y) (see ( Corresponding theorems to Theorems 3.1, 3.2, 3.3 and 3.5 then follow with some changes to their sufficient conditions. Analogous comparison results to Theorems 4.1 and 4.2 hold. Further, if all mixing events are of size c and π c = min(ζ /c, 1), where ζ ∈ (0, ∞), then if R 0 > 1 and γ are held fixed and c → ∞, the probability of a minor outbreak given one initial infective, z, tends to a limit that is strictly less than 1, while the size of a major outbreak, τ , tends to 1 − R −1 0 . The limit of τ for the corresponding model with Reed-Frost type mixing is strictly greater than 1 − R −1 0 . A few models with non-pairwise transmission are mentioned in the introduction. Other examples include models incorporating synergystic interactions (eg. Ludlam et al. 2012), models in which the usual mass-action law of infection, where the overall force of infection is proportional to the product of the numbers of susceptible (S) and infective (I ) infective individuals, is replaced by some other function of (S, I ) (eg. O'Neill and Wen 2012 and references therein) and in nonparametric inference, such as Knock and Kypraios (2014), where no particular form is assumed for the force of infection. In all of these models new infectives occur one at a time. Billard et al. (1980) analyse an SI model (i.e. one in with no recovery from infection) in which infections occur in batches whose sizes are i.i.d. random variables. However, models in which infection is not driven by pairwise interaction of individuals are relatively rare.
The model was motivated by non-pharmaceutical intervention policies for Covid-19 that placed limits on the size of gatherings outside the home and it would be interesting to use it to explore the effects of such policies. It is straightforward to analyse, at least numerically, the effects of changes in the distribution of the mixing group size random variable C on epidemic properties such as R 0 and the probability and size of a major outbreak. The LLN and functional CLT in Sect. 3.2 can be extended to allow the distribution of C to be time and/or state dependent, thus permit-ting investigation of corresponding control measures. Another feature that is highly pertinent to Covid-19 spread and control is the role played by households. It would be worthwhile to examine the effects of extending the households model of Ball et al. (1997) so that between-household spread is modelled using mixing groups, rather than by homogeneous mixing. Another promising direction for future work is to introduce heterogeneities, eg. owing to age and/or level of social activity, and study the corresponding multitype model.