Herd immunity under individual variation and reinfection

We study a susceptible-exposed-infected-recovered (SEIR) model considered by Aguas et al. (In: Herd immunity thresholds for SARS-CoV-2 estimated from unfolding epidemics, 2021), Gomes et al. (In: J Theor Biol. 540:111063, 2022) where individuals are assumed to differ in their susceptibility or exposure to infection. Under this heterogeneity assumption, epidemic growth is effectively suppressed when the percentage of the population having acquired immunity surpasses a critical level - the herd immunity threshold - that is lower than in homogeneous populations. We derive explicit formulas to calculate herd immunity thresholds and stable configurations, especially when susceptibility or exposure are gamma distributed, and explore extensions of the model.


Introduction
Understanding and predicting the dynamics and control of infectious diseases relies on representative models, whether conceptual or mathematical. Mathematical modelling was established in infectious diseases over a century ago, with the seminal works of Ross (1916), Ross and Hudson (1917), Kermack and McKendrick (1927) and others. Propelled by the discovery of aetiological agents for infectious diseases, and the germ theory, models have focused on the complexities of pathogen transmission and evolution (Heesterbeek 2015). It has recurrently been noted for over a century, however, that these models tend to overpredict transmission potential and overestimate the impact of control measures which may be explained by limitations in capturing the effects of heterogeneity (Kermack and McKendrick 1927;McKendrick 1939;Gart 1968Gart , 1971Ball 1985;Anderson et al. 1986;Pastor-Satorras and Vespignani 2001;Miller et al. 2012;Gomes et al. 2022).
Here we analyze a set of susceptible-exposed-infected-recovered (SEIR) models presented in Gomes et al. (2022), Aguas et al. (2021) where each of the compartments S, E, I and R is expanded into continuum many compartments S(x), E(x), I (x) and R(x), where x ∈ R + is a trait that varies among individuals. Specifically we model a situation where each individual has a level of susceptibility or exposure (connectivity) x, starting in compartment S(x) and staying within the compartments S(x), E(x), I (x) and R(x) the whole time. This individual may infect or be infected by others irrespective of their trait value x assuming random mixing (Anderson and May 1991;Diekmann et al. 2013). We will consider two types of models: A variable susceptibility case where the susceptibility of an individual at level x is proportional to x, or, in other words, if we compare an individual at level x and an individual at level y, the one at level x is x/y times more likely to get infected than the one with susceptibility y. We may interpret this as variation in biological susceptibility which may be due to genetics, epigenetics or life history.
A variable connectivity case where the propensities for an individual at level x to acquire infection and transmit to others are both proportional to x, or, in other words, if we compare an individual at level x and an individual at level y, the one at level x is x/y times more likely to get infected than the one in level y and also x/y times more likely to infect someone else once infected. This is interpreted as individuals with many contacts being both more likely to get infected and to infect others.
For each x, we have a system of the form: where λ is the force of infection which is formulated differently in the variable susceptibility or the variable connectivity cases: Variable susceptibility: λ = β I (x) dx, Variable connectivity: λ = β x I(x) dx.
Note that λ varies with time, as it depends on the time-dependent infected population. The dynamics of compartments S(x), E(x), I (x) and R(x) are governed by the infinite system of ordinary differential equations: We assume that the system has been scaled such that the total population is 1. The initial conditions for variables S(x, t), E(x, t), I (x, t) and R(x, t) 1 is a small scalar to seed the epidemic, and q(x) is a probability density function with mean 1 and coefficient of variation ν: We use S(t) to denote the integral over all susceptibility levels, S(x, t), for x ∈ R + . We thus have S(t) = +∞ 0 S(x, t)dx. Same with E(t), I(t) and R(t).
We will use the first three moments of S(x, t), that we denote S(t),S(t) and S(t): When infection is absent ( = 0), we have S(0) = 1,S(0) = 1 and S(0) = 1 + ν 2 . But note that S(x, t) is not a probability density function for > 0 as S(t) becomes less than 1. The quotient S(x, t)/S(t) as a function of x for fixed t will be a probability density function for > 0 and all t with first and second momentsS(t)/S(t) and S(t)/S(t) which decrease over time. When the initial configuration q(x) is a gamma distribution, all the distributions S(x, t)/S(t) are also gamma with the same coefficient of variation ν but with lower mean (see Appendix A), an argument which enables mathematical derivations to advance further when traits are assumed to be gamma distributed (Novozhilov 2008).
Similarly, we define the moments R(t),R(t) and R(t) for the recovered compartment, and the same with E and I. Notice for instance that λ(t) is written as β I(t) and βĪ(t) in the variable susceptibility and variable connectivity cases, respectively.
Here we describe key epidemiological quantities when system of equations (Eqs. 1-4) is adopted. The basic reproduction number R 0 is the average number of secondary infections generated by an infected individual in a totally susceptible population. It depends on characteristics of both the pathogen and the host population. When this number is below 1 no epidemics are expected. When R 0 is above 1, however, the introduction of infection in a virgin population is expected to generate an epidemic. This is followed by almost exponential growth in cumulative infections which decelerates gradually as susceptibles are depleted. The effective reproduction number R eff is a time-dependent quantity loosely defined as the number of secondary infections generated by a typical infected individual when the susceptibility of the population is as at time t. R eff coincides with R 0 at the beginning of an epidemic (when the population is totally susceptible) but declines as individuals are removed from the susceptible pool by infection and immunity. As R eff crosses 1 towards lower values, the epidemic subsides and future reintroductions of infection are not expected to generate new outbreaks as long as population immunity is maintained.
We derive formulas for the effective reproduction number R eff and the herd immunity threshold H in terms of moments S,S and S, of the susceptible population. When q(x) is a gamma distribution,S and S can be formulated in terms of S and we get an exact formula for H in terms only of the basic reproduction number R 0 and the coefficient of variation ν. In this case we can also reduce the infinite system (Eqs. 1-4) to a finite system of ordinary differential equations in S, E, I and R with nonlinear transmission (exactly when the variable trait is susceptibility and approximately in the case of variable connectivity). In the case of variable connectivity, we provide an exact derivation of a finite system in the variablesS,Ē,Ī andR.
In the variable susceptibility case we will get that: and consequently, R eff = R 0S . This implies that the population is above the herd immunity threshold whenS < 1/R 0 . When we assume that q(x) is a gamma distribution, the proportion of individuals that have been infected by the time the herd immunity threshold is reached is deduced as: and system (Eqs. 1-4) can be reduced to: These equations are exact. In the variable connectivity case we will get that: and consequently, R eff = R 0 /(1 + ν 2 ) S. This implies that the population is above the herd immunity threshold when S < (1 + ν 2 )/R 0 . When we assume that q(x) is a gamma distribution, the proportion of individuals that have been infected by the time the herd immunity threshold is reached is deduced as: In (Eqs. 28-31), we derive a closed system for the variable connectivity model on the variablesS,Ē,Ī andR. But this does not directly allow the model to be fitted to incidence data provided by routine surveillance. For finding equations that determine the system only on the variables S, E, I and R, in the variable connectivity case we need to make an approximation as detailed in Appendix B. The resulting system (Eqs. 22-25) is shown to approximate the original (Eqs. 1-4) when the infectious period is small as in acute infectious diseases. In Fig. 1 we provide graphical representations for the H formulas corresponding to the basic model introduced in this section (variable susceptibility in green and variable connectivity in blue), showing a monotonic decrease as the coefficient of variation ν increases (Gomes et al. 2022).
In Sects. 2 and 3 we provide general derivations for effective reproduction numbers and herd immunity thresholds, while Sect. 4 is focussed on special cases when traits x are gamma distributed. Towards the end of the paper, we analyze two extensions of the basic model. In Sect. 5, we consider a model with reinfection where immunity after recovery is not fully protective but only partially. In Sect. 6, we consider a model with a carrier state, which (Gomes et al. 2022;Aguas et al. 2021) apply to the coronavirus disease (COVID-19) pandemic. There, the exposed compartments are not simply latent but a carrier state where individuals are infectious but to a lesser degree than individuals in the fully infectious compartment. In each case we derive formulas for herd immunity thresholds especially when the initial trait distribution is gamma.

Effective reproduction number
The effective reproduction number at time t is defined as the number of secondary infections caused by a typical infected individual over their entire infectious period in an idealized situation, as described below.
Before proceeding with its derivation we need to make two considerations. The first concerns the evolving susceptible pool. We assume that during the period in which an individual is contagious, the density of susceptibles is frozen in time. This means that we disregard the fact that since the susceptible population declines, this individual infects less at the end than at the beginning of their infection. In acute infections, the decline in the susceptible population is usually slow compared to the rate of recovery from infection so, in practice, the impact of the assumption is negligible. Moreover, when R eff is used to analyze stable configurations, such as in the derivation of herd immunity thresholds, the assumption holds and hence has no effect on the results. This consideration pertains to both variable susceptibility and variable connectivity models.
The second concerns the infectivity profile of the infected population at time t. We define: R eff at time t as the average number of secondary infections generated by an individual who becomes infected at time t. This average is taken over the pool of individuals that go from S to E at time t.
When R eff < 1, infection is not expected to invade an infection-free population. Further details on this concept are discussed in Appendix B. In the remaining of this section we derive explicit formulas for R eff .
First, the variable susceptibility case: Consider an individual who gets infected (more precisely, exposed and consequently infected) at time t (i.e., moves from S to E at time t). This individual will eventually move to I and spend on average 1/γ days there. While in I, the individual will infect an average of β y S(y, t) dy others per day. We thus get: In particular, we get: and consequently: Second, the variable connectivity case: Consider again an individual who gets infected at time t. It now matters what trait value x this individual has because it determines how many others they will infect.
Let p(x, t) be the density function measuring the probability at time t that this individual has connectivity level x. The probability of becoming infected (i.e., of As above, an individual who enters E will eventually move to I, spend on average 1/γ days there, and infect an average of β y S(y, t) dy others per day. We thus get: In particular, we get: and consequently: Expressions for R 0 , such as (Eq. 14) and (Eq. 17), have been known for decades (Anderson and May 1991;Diekmann et al. 2013;Woolhouse et al. 1997). Worth highlighting, however, is that disproportionately less attention has been given to variable susceptibility than to variable connectivity due to the coefficient of variation ν not affecting the formula explicitly in the former case but only in the latter. It will be important to realise, however, that variation in susceptibility is just as impactful when we consider quantities such as herd immunity thresholds and inferences of R 0 from observational data (Gomes et al. 2022;Aguas et al. 2021).

Herd immunity threshold
Suppose we have a population with no infected individuals, so that all individuals are either susceptible or recovered. The population is said to be at or above the herd immunity threshold for a pathogen if its susceptibility profile to that pathogen is such that a new introduction of infection (i.e., a small increase in E or I) does not trigger an outbreak. By inspection on the differential equations (Eqs. 1-4), we see that a configuration with E(x) = I (x) = 0 satisfies this condition if and only if R eff ≤ 1.
In the variable susceptibility case it is equivalent to formulate the herd immunity threshold in terms of suppression of future outbreaks (as adopted here) or in terms of an unmitigated epidemic passing its peak infection prevalence. With variable connectivity, however, this equivalence does not hold as explained in Sect. 2.
A configuration with no infected individuals is then said to be at the herd immunity threshold if and only if R eff = 1. In SEIR models with no individual variation, configurations with no infected individuals are determined by the values of S and R = 1 − S, and the herd immunity threshold is defined as the value of 1 − S at the unique configuration with R eff = 1. This value is well-known to be equal to 1 − 1/R 0 (Anderson and May 1991;Diekmann et al. 2013). With individual variation in susceptibility or exposure to infection, however, there are many configurations which satisfy R eff = 1. One such configuration is given by S(x) = q(x)/R 0 for all x. This could be obtained, for instance, by vaccinating a proportion 1 − 1/R 0 of the total population randomly without taking into account susceptibility or exposure levels (Fine et al. 2011). When immunity is acquired naturally, however, individuals with higher susceptibility or exposure tend to be infected earlier and the herd immunity threshold is reached before the susceptible population is as low as 1/R 0 of the total. For now, let us retain that for the basic heterogeneous models considered here the herd immunity threshold is reached whenS = 1/R 0 in the variable susceptibility case (Eq. 15) and S/(1 + ν 2 ) = 1/R 0 with variable connectivity (Eq. 18).
Next, we will see how we can deriveS(t) and S(t) under the assumption that the initial distribution is gamma.

Case of the gamma distribution
Here we study how the distribution of the trait x within the susceptible compartment evolves generally, and then specify to the case where individual variation in susceptibility or connectivity is gamma distributed. We also refer to related work in Neipel et al. (2020), Novozhilov (2008.

Evolution of the susceptible compartment
From the SEIR equation for d S(x, t)/dt we get: Integrating with respect to t we get: This holds in both variable susceptibility and variable connectivity models (with different values for k t ). It also holds in the cases with reinfection and with carrier state considered in Sects. 5 and 6.

Gamma distributed traits
The key observation here is that since S(x, t) = q(x) e −x k t , we have that S(x, t)/S(t) remains a gamma distribution at all values of t. This enables the derivations in Appendix A of explicit formulas for the momentsS (Eq. A.1) and S (Eq. A.2) in terms of S and the shape parameter α, when susceptibility or connectivity are gamma distributed.
Using that the coefficient of variation is ν = 1/ √ α we rewrite the respective formulas as:S To derive the reduced system in the variable susceptibility case (Eqs. 8-11) we integrate the equations in (Eqs. 1-4) over the susceptibility domain and apply (Eqs. 13 and 20) to get: This closed system in S, E, I and R has been used to fit epidemic curves of COVID-19 (Gomes et al. 2022;Aguas et al. 2021). Recalling that the herd immunity threshold H is 1 − S(t) at the time t when R eff = 1 and that R eff (t) = R 0S (t) (Eq. 15), we get that R −1 0 =S(t) = S(t) 1+ν 2 at the herd immunity threshold when susceptibility is gamma distributed. Hence: In the variable connectivity case we have that R eff (t) = R 0 /(1+ν 2 ) S(t) (Eq. 18). Thus, when R eff = 1, we get R −1 0 = S(t)/(1 + ν 2 ) = S(t) 1+2ν 2 when connectivity is gamma distributed. Hence: Multiplying the original (Eqs. 1-4) by x, integrating over the connectivity domain and applying (Eqs. 16 and 21) we get: Mathematically this is a tractable closed system inS,Ē,Ī andR. However, these variables are not convenient for practical data fitting and parameter estimation. In Appendix B we propose an approximation in the variables S, E, I and R.

Model with reinfection
Here we consider an extension of the model considering that immunity after recovery is not fully protective, but only partially. A factor σ , with 0 ≤ σ ≤ 1, is added to represent the quotient of the probability of getting reinfected after recovery over the probability of getting infected while fully susceptible.
The model is represented diagrammatically as: with λ as in the basic models (without reinfection) studied above. The extended model is given by the equations: The system exhibits newer dynamics in comparison with the basic case. Depending on whether σ is below or above 1/R 0 (known as the reinfection threshold (Gomes et al. 2004(Gomes et al. , 2016) we get that either the disease dies out after a while and a certain proportion of the population never gets infected, or continues endemically and every individual is eventually infected and reinfected repeatedly.

Effective reproduction number
The basic reproduction number is calculated exactly, as in the absence of reinfection, but the effective reproduction now depends not only on the distribution of S(x, t), but also on the distribution of R(x, t). When we consider configurations with no infected individuals, we will have that R(x, t) = q(x) − S(x, t) and will be able to express R eff (t) in terms of S(x, t) only.
The derivation of these formulas is essentially as the derivations in Sect. 2, with two differences. First, each individual with trait value x infects β(S(t) + σR(t)) or βx(S(t) + σR(t)) others per day spent in I, in the respective cases, instead of βS(t) or βxS(t). Second, when we consider an individual that gets infected in the variable connectivity case, the probability that this individual has trait value x is proportional to x(S(x, t) + σ R(x, t)) instead of x S(x, t).

Herd immunity threshold
Recall that a configuration with no infected individuals is at or above the herd immunity threshold if and only if R eff ≤ 1. Assuming that no one is infected, that is R(x) = q(x) − S(x), we getR = 1 −S and R = q − S, where q = x 2 q(x) dx = 1 + ν 2 . We can then understand the configurations at the herd immunity threshold in terms of S and S.
In the variable susceptibility case a configuration with no infected individuals is at the herd immunity threshold if and only ifS + σ (1 −S) = 1/R 0 , and hence: In the variable connectivity case a configuration with no infected individuals is at the herd immunity threshold if and only if (S + σ (q − S))/(1 + ν 2 ) = 1/R 0 , and hence:

Reinfection threshold
The formulas above require: That is, the reinfection factor σ has to be below R −1 0 , a critical value known as the reinfection threshold (Gomes et al. 2004(Gomes et al. , 2016. If this is verified, then all configurations with no infected individuals and satisfying the conditions above (either depending on the case) are herd immunity threshold configurations in the sense that when susceptibility is at that level or lower, infection reintroductions will not trigger new outbreaks as R eff will not increase above 1.
If the reinfection factor is so high that the population is above the reinfection threshold, R eff will be greater than 1 in any such configurations, so there will not be any configuration with no infected individuals which is at the herd immunity threshold. This implies that there will always be a portion of the population infected, and hence that the population of susceptible individuals will eventually be completely depleted. The infection becomes endemic. The equilibrium configuration will now have S(x) = 0 for all x. In these situations the level of endemicity will depend on how much resistance the population is able to mount and maintain.

Case of the gamma distribution
Recall from Sect. 3 that, when the initial distribution q(x) is gamma, we get that S(x, t)/S(t) remains a gamma distribution for all t, and thatS(t) = S(t) 1+ν 2 and S(t) = (1 + ν 2 ) S(t) 1+2ν 2 . We can then obtain the values of S(t) at the time when the herd immunity threshold is reached, and then calculate H as 1−S(t) for that particular t.
In the variable susceptibility case we get: In the variable connectivity case we get: Curves assuming a selection of values for σ are represented graphically in Fig. 2. Note the critical behaviour at the reinfection threshold (σ = 1/R 0 ) in red, which separates the regime where individual immunity is sufficiently potent for the herd immunity threshold to be achievable from the regime where endemicity will establish. Herd immunity threshold with reinfection. Curves correspond to the SEIR model with reinfection (Eqs. 32-35) assuming R 0 = 3 and gamma distributed susceptibility (left, Eq. 36) or connectivity (right, Eq. 37). Efficacy of naturally acquired immunity is captured by a reinfection parameter σ , potentially ranging between σ = 0 (100% efficacy) and σ = 1 (0 efficacy). Five values of the reinfection parameter are depicted: σ = 0 (black); σ = 0.1 (green); σ = 0.2 (blue); σ = 0.3 (magenta); and σ = 1/R 0 (red). Above σ = 1/R 0 (reinfection threshold Gomes et al. 2004Gomes et al. , 2016 the infection becomes stably endemic and there is no herd immunity threshold (colour figure online)

Model with a carrier state
In applications to the COVID-19 pandemic (Gomes et al. 2022;Aguas et al. 2021) the exposed compartments are not simply a latent state but a carrier state where individuals are infectious but to a lesser degree than those in the fully infectious compartment. Building on the reinfection model (Eqs. 32-32) we now introduce parameter ρ ≤ 1 to denote the ratio of infectiousness between exposed and fully infectious individuals. What changes is the force of infection λ: Variable connectivity: λ(t) = β x(ρ E(x, t) + I (x, t)) dx = β (ρĒ(t) +Ī(t)).

The basic and effective reproduction numbers become:
Variable susceptibility: Variable connectivity: The calculation of the effective reproduction number R eff is slightly different. The difference is that now we have to add the time an individual is incubating the infection in E to the infectious period, multiplied by the factor ρ. Since the average time an individual spends in E is 1/δ we get: • R eff (t) = β(ρ/δ + 1/γ ) (S(t) + σR(t)) in the variable susceptibility case, and • R eff (t) = β(ρ/δ + 1/γ ) (S(t) + σ R(t)) in the variable connectivity case, whereS(t) and S(t) are the moments of S(x, t) defined above, and the same with R(x, t).

Discussion
After completion of this work, similar ideas for capturing individual variation with mean-field epidemic models were elaborated (Britton et al. 2020;Di Lauro et al. 2021;Kawagoe et al. 2021;Neipel et al. 2020;Rose et al. 2021;Tkachenko et al. 2021). These recent developments were largely prompted by the COVID-19 pandemic. In agreement with our results, Rose et al. (2021), Neipel et al. (2020) find that when susceptibility is initially gamma distributed it remains so through the course of the epidemic, leading naturally to power-law behaviour in the force of infection (Novozhilov 2008). In addition the authors show that other initial distributions converge towards gamma through the process of contagion. Other authors (Kawagoe et al. 2021) derive epidemic final sizes assuming alternative distributions of susceptibility, considering in addition that infectivity may exhibit some correlation with susceptibility (such as in the variable connectivity models analyzed here). They compare numerical results for gamma and lognormal distributions with those obtained using an empirical distribution of individual contacts derived from cell phone geolocation data. Focussing on variable connectivity (Britton et al. 2020;Di Lauro et al. 2021) address herd immunity thresholds using age-structured compartmental models. Additionally (Di Lauro et al. 2021) consider a variety of non-pharmaceutical intervention scenarios, emphasizing subtle results when interventions change the contact network. Another group of authors (Tkachenko et al. 2021) distinguishes between persistent and transient individual variation to highlight that only the former is subject to the kind of selection that lowers epidemic final sizes and herd immunity thresholds.
Despite being prompted by COVID-19 none of the above references attempted to quantify the individual variation which was under selection by the force of infection and hence contributed to lower epidemic final sizes and herd immunity thresholds. This is done in our associated work (Gomes et al. 2022;Aguas et al. 2021), where the problem is inverted and selectable variation is inferred from its effects on epidemic patterns. Although we and others had previous adopted similar approaches to other infectious diseases (Finkenstadt and Grenfell 2000;Smith et al. 2005;Bellan et al. 2015;Corder et al. 2020) the work has been processed more cautiously during the pandemic due to the greater implications that estimating lower herd immunity thresholds and epidemic final sizes might have for policies and behaviors in this context. The concept of herd immunity was originally developed in the context of vaccination programs (Gonçalves 2008;Fine et al. 2011). Defining the percentage of the population that must be immunised to cause infection prevalence to decline, the concept has provided a useful target for vaccination coverage. In hypothetical scenarios of vaccines delivered at random and individuals mixing at random, the herd immunity threshold is given by the simple expression (1 − R −1 0 ) when immunity is fully protective. Concretely, for R 0 between 2.5 and 5, this would indicate that 60-80% random subjects would need be immunized to prevent spread of infection. This formula would not apply, however, if vaccination programmes were designed to prioritize more connected individuals, for instance (Elbasha and Gumel 2021). Similarly, it does not apply when immunity is acquired in response to natural infection, which does not occur at random. Individuals who are more susceptible or more exposed to infection are prone to be infected and become immune earlier. As a result, earlier episodes contribute disproportionately to herd immunity as they remove highly susceptible and exposed subjects from the susceptible pool (Ferrari et al. 2006;Britton et al. 2020;Kawagoe et al. 2021;Neipel et al. 2020;Rose et al. 2021;Gomes et al. 2022;Aguas et al. 2021).
In our basic models, the herd immunity threshold becomes H = 1 − R −1/(1+ν 2 ) 0 in the case of gamma distributed susceptibility, and H = 1 − R −1/(1+2ν 2 ) 0 with gamma distributed connectivity (exposure), which decline sharply when coefficients of variation (ν) increase from 0 to 2, remaining below 20% for more variable populations in a particular illustration where R 0 = 3 (Fig. 1). The magnitude of the decline depends on what property is heterogeneous and how it is distributed among individuals, but the downward trend is robust provided that acquired immunity is efficacious enough to keep transmission below the reinfection threshold ( Fig. 2) (σ < R −1 0 , where σ is the susceptibility of individuals who have recovered relative to their respective susceptibility prior to infection) (Gomes et al. 2004(Gomes et al. , 2016. In our reinfection models, herd immunity thresholds are derived as in the case of gamma distributed susceptibility, and with gamma distributed connectivity, when σ < R −1 0 . If immunity is not potent enough to keep the system below the reinfection threshold then a herd immunity threshold is not attainable and the disease persists in stable endemicity, irrespective of individual variation. Finally, we stress that the herd immunity threshold (H) is a theoretical framework to assess epidemic potential to the same extent that R 0 is a theoretical framework. Their interdependence shows that if R 0 increases due to evolution of the infectious agent, for example, so does H. Also, if new susceptibles enter the population through birth or other processes, or if immunity wanes or is evaded by pathogen lineages, a previously acquired herd immunity status may be lost leaving the population prone to subsequent outbreaks. Furthermore, if transmission has a marked seasonal pattern the same level of immunity may place the population above threshold in low season and below threshold in high season, in a cyclical manner. Although H is not as immediately applicable as often implied, it is a more informative measure of epidemic potential than R 0 given that it accounts for variation in susceptibility or exposure (ν) in addition to average transmissibility R 0 . The more accurately we know H the better we can assess trade-offs and inform public health policy. and S(t): From the above we getS

Appendix B. Approximate variable connectivity model in S, E, I, R variables
In (Eqs. 28-31), we derived a closed system for the variable connectivity model on the variablesS,Ē,Ī,R. However, we would like to have a closed system on the variables S, E, I, R, as we did in the variable susceptibility case to enable direct fitting to incidence data  for parameter estimation and scenario projections. We propose a system of approximate equations which works well when the infectious period is small as in acute infectious diseases (Eqs. B.8-B.11 below). The justification for the approximation is based on the assumption that the following two quantities are very close to each other when the infectious period is small compared to the length of the epidemic. At each time t we consider the following two factors, that we denote f(t) and t(t) • f(t) is the average infectivity of a typical individual who becomes infected at time t.
In other words, this is the average value of trait x taken over the pool of individuals that go from S to E at time t. • t(t) is the average infectivity of a typical individual who is infectious at time t. In other words, this is the average value of trait x taken over the pool of individuals that are in compartment I at time t.
Both f and t are formulated below (Eq. B.2). The distinction between them is only meaningful in models where infectivity varies among individuals in such a way that those infected in different days tend to have different infectivities -this occurs in our variable connectivity model (as infectivity is positively correlated with susceptibility, the selection which reduces susceptibility in the susceptible pool over time also reduces infectivity) but not in the case of variable susceptibility (where all individuals have the same infectivity). Note that we use the factor f in our calculation of R eff in (Eq. 16), where we multiply β by the average infectivity f (same as average connectivity), the average susceptibilitȳ S, and the length of the infectious period (1/γ ): Let us also note that the "R-number" which was popularized in the COVID-19 pandemic, loosely described as the number of new infections caused by a typical individual who is infectious at time t, is defined in terms of t rather than f. This more empirical notion, formulated as R t = (β/γ ) · t ·S, is not always coincident with the basic reproduction number R eff , introduced in textbooks as the spectral radius of the next generation operator (Diekmann et al. 2013) and adopted in this paper. The two notions differ when the average infectivity of individuals infected in a particular day is different to the average infectivity of those who have been infected in another day. However, this difference is of little consequence to variable connectivity models when the infectious period is small and of no consequence when variation is in susceptibility (see Fig. 3). We claim that the values of f and t are given by the following formulas: The second equality follows from the formula for average infectivity, (x I (x, t)/I)dx, for every time t. The first equality was essentially derived in Sect. 2 when we calculated R eff for variable connectivity (Eq. 16). Within that calculation we had to consider the infectivity of an individual who becomes infected at time t, which is equal to x p(x, t)dx, where p(x, t) is the density function measuring the probability that such an individual has connectivity level x. We saw that since p(x, t) is proportional to x S(x, t), we have that p(x, t) = x S(x, t)/S(t). It follows then that f(t) = x 2 S(x, t)/S(t)dx, for every t, which is the expression above. We now can write our assumption as: With this we can derive the approximate closed system of equations on the variables S, E, I, R. Observe that if we integrate the Eqs. Top panels show epidemic curves generated by running system (Eqs. 1-4) with explicit gamma distributed susceptibility (left) or connectivity (right). Solid curves are unmodified epidemics while dashed show the outcome of hypothetically clearing all infections instantaneously on day 25 and progressing with a small infectious seed. Bottom panels show the respective R t values, calculated as (β/γ ) ·S in the case of variable susceptibility and (β/γ ) · (Ī/I) ·S in the case of variable connectivity, with or without the intervention on day 25. The empirical index R t coincides with the effective reproduction number R eff in the variable susceptibility case but is slightly higher when variation is in connectivity. This is highlighted by zooming into the areas marked by black rectangles (see insets). Parameters: R 0 = 3; δ = 1/4 per day; γ = 1/3 per day; and ν = 1 (colour figure online) Then, using thatĪS ≈ I S (Eq. B.3), and that with gamma distributed traits we have S = (1 + ν 2 ) S(t) 1+2ν 2 (Eq. 21), we get the approximate closed system we seek: illustrates epidemic curves generated with system (Eqs. 1-4) with an explicit gamma distribution (mean 1 and coefficient of variation ν = 1) for susceptibility (solid green) and connectivity (solid blue). Technically the distribution was discretized in n bins and a system of 4n ODEs was run to generate the curves. Dashed curves were added to the plots to enable comparisons with the corresponding results generated by the reduced systems respectively). The bottom panel shown the associated effective reproduction numbers R eff in the case of the reduced systems, and the empirical proxy R t when explicit distributions are implemented (recall that the two coincide when variation is in susceptibility). All runs were conducted with basic reproduction number R 0 = 3, incubation period (1/δ = 4 days), infectious period (1/γ = 3 days) and coefficient of individual variation (ν = 1). Notice that the reduction is exact in the variable susceptibility case but not with variable connectivity. In a related study we fitted both explicit and reduced models to real epidemics and obtained similar parameter estimates and conclusions irrespective of which version was adopted.