The recent emergence of SARS-CoV-2 and the associated disease COVID-19 has had worldwide impact. Many cities have had large outbreaks and brought them under control through major interventions. Once those interventions are lifted, in absence of effective vaccination, a homogeneous model of infection spread would predict that as long as less than 1 − 1/R 0 of the population was infected, there is always a threat of resurgence.
Despite large epidemics, cities such as New York remain well below the threshold expected to be be required to achieve herd immunity [27] .
To avoid the significant economic and health costs associated with continued interventions, it is natural to consider the so-called "herd immunity" strategy. This strategy allows infection to spread with restrictions in place so that the outbreak finishes and interventions are lifted when the herd immunity threshold is reached. Typically in an uncontrolled epidemic the herd immunity threshold is reached at the peak of the epidemic, and many additional infections occur as the outbreak slowly wanes. The additional infections are sometimes referred to as the overshoot [11] . By calibrating the intervention so that there are no (or almost no) infections when the herd immunity threshold is reached, interventions can be removed with minimal overshoot [7, 20] .
The severity of the epidemic in many places whose seroprevalence is still very low has led many to suggest that the herd immunity strategy is not tenable. However, recent papers [3, 10] suggest that immunity acquired through infection may be distributed through the population in such a way as to achieve herd immunity at a lower fraction affected than a homogeneous model would predict. This is because the initial wave of infections preferentially affects those most at risk. Thus it acts like a targeted vaccination, removing the people who are most likely to transmit infection from the susceptible pool. Generally, the time course of a real epidemic involves the following stages: (i) short period of unconstrained transmission, (ii) significant control or lockdown, and (iii) relaxation of control measures. Typically, during lockdown some spread persists and one pertinent question is whether relaxing the lockdown will lead to a second wave. In [3] , this question was explored by looking at prevalence at the end of lockdown in scenarios where lifting the lockdown did not lead to resurgence. For the purposes of this paper, we refer to this fraction as the disease-induced herd immunity (DIHI).
These papers make some simplifying assumptions about population structure that may not hold. In particular they do not consider the fact that existing interventions tend to affect some contacts more than others. For example, the transmission rates of household contacts are not significantly reduced (and in fact may increase) during interventions focused on reducing movement. Moreover, many of the highest risk positions in the disease network are in fact roles (health care workers, delivery drivers, teachers, etc) that would need to be maintained in most forms of a lockdown, meaning that some community links may be increased even as others are limited.
In this paper we will use several network-based models to improve our understanding of how population structure may affect the DIHI threshold. Modelling epidemics on networks and the analysis of the resulting systems, be they stochastic or deterministic, is a well studied research area [14, 23, 25] . The main reason for using networks is twofold. First, networks are intuitive and easy to understand and use, and second, they allow us to represent contact structure within a population to a high degree of accuracy. One of the most frequently studied questions in this area is how to understand and quantify the impact of various network properties on the invasion, spread and control of infectious diseases. Epidemics are usually modelled as stochastic process unfolding on the network. For example, the susceptible-infected-recovered/removed (SIR) epidemic on a network is modelled by the two separate processes: infection and recovery. In the simplest case both are Poisson processes such that an infectious node u connected to a susceptible node v transmits at per-contact rate τ , causing v to convert to infected. Nodes recover at rate γ to an immune state independently of the network.
On a network with N nodes, the SIR model leads to a continuous-time Markov Chain over 3 N states. Even for small values of N such a system becomes tedious to handle (numerically or otherwise). This has led to a myriad of approximations in the form of mean-field models which focus on some average quantities to reduce the dimensionality of the resulting system [14] .
The mathematical analysis of epidemic models on networks (from regular and Erdős-Rényi to scalefree networks such as the Barabási-Albert preferential attachment model) has led to many elegant analytical / explicit results regarding the impact of network properties (e.g., degree heterogeneity, assortativity and clustering) on the epidemic threshold, final epidemic size or endemic equilibrium, optimal vaccination strategies [25, 23, 5, 12] , as well as the fact that certain mean-field SIR models are exact in the limit of large configuration model networks [19, 17] . Many such models are mathematically tractable and help us gain important intuition about how an epidemic spreads and what properties of the contact network affect this. While they may not directly inform policy making, they remain useful tools to develop intuition, highlight important population features and test various scenarios quickly and effectively.
In this paper we use four different mean-field models [14] to approximate exact epidemics on networks: (a) degree-based heterogeneous mean-field, (b) clustered and unclustered pairwise, (c) a new edge-based-compartmental model that allows us to distinguish between household/local contacts and community/global contacts, and (d) an age-structured compartmental model parameterised with realistic age structure and contact matrices. The primary aim is to investigate how heterogeneity in model structure impacts DIHI. More importantly perhaps, we challenge the way lockdown has been implemented in many models, namely, by a simple reduction in R 0 or the transmission rate while keeping the contact network or mixing matrix the same. We build a new edge-basedcompartmental model able to distinguish between household and community transmissions and use this to implement lockdown by either intervening on both types of connections or only on the community-based ones. In the same spirit, we use an age-structured compartmental model in which we implement lockdown either as a simple scaling of the entire mixing matrix or, more realistically, a set of age-specific structural changes.
The paper is structured as follows. In section 2 we describe all models including the underlying network types and relevant model and epidemic parameters. Section 3 contains the results for the network models, whilst section 4 provides the results for the age-structured model. Finally, in section 5 we discuss the implications of our findings. Additional technical details are given in the Appendix.
We consider a set of mean-field models. The first two are approximations of the exact stochastic SIR model on networks with heterogeneous degree distribution, without and with clustering. This is followed by a new edge-based-compartmental model with household structure, i.e., with the ability to distinguish between household and community contacts. Finally, we consider an age-structured SEIRD model based on realistic age-structure and mixing matrices. For the network-based meanfield models we choose a flexible degree distribution with good control of the mean and variance, i.e., the negative binomial. We loosely order the models by their relative complexity, corresponding to gradually incorporating more features of the underlying network or population contact structure. While the last model is not explicitly network-based, it uses realistic mixing matrices over an age-structured model based on UK data. One can consider this an extrapolation from (or an approximation of) an explicit network model where individual-level interactions are averaged out over groups of interest from an epidemiological view point.
We consider a SIR/SEIRD epidemic spreading in a closed population of size N = 6.65 · 10 6 (loosely the population size of the UK) with a well-defined contact structure. For illustrative purposes, we assume that the probability of an individual having k contacts follows a negative binomial distribution
The reason for this choice is that we want to highlight how heterogeneities in the contact structure play a central role in determining the DIHI. To illustrate this point, we consider three different scenarios for the degree distribution of the population. In all cases, we fix k = n(1 − p)/p and we use the remaining free parameter to tune the variance. To avoid individuals with degree 0, the degree distribution is shifted by a constant m, thus making the effective average degree k = m + n(1 − p)/p. For normal-like and scale-free-like distributions we take m = 1, and for the delta-like distribution we take m = 9 (see Table 1 ).
The parameters chosen are reported in Table 1 , and represent degree distributions of increasing variance, see figure 1 , moving from almost no variance (delta-like degree distribution) to a degree distribution with a longer tail, akin to a scale-free network.
For all network-based models, we consider the simple susceptible-infected-recovered model (SIR) on a network. The infection is driven by a per-link transmission rate τ and a recovery rate γ. Mostly, we focus on the corresponding mean-field models with increasing level of accuracy in terms of what network features are incorporated. Table 1 : The three degree distributions considered. The delta-like distribution is shifted by 9, as its mean would be 1 otherwise. The reason for this choice is the fact that in negative binomials the variance cannot be lower than the mean. Normal-like and scale-free like distributions instead are shifted by one, so that the minimum degree is 1. The resulting degree distributions are shown in figure 1 . Table 1. meaning that the infection pressure across a link is simply averaged across the entire spectrum of infected nodes. The resulting ODEs arė
where N = P n,p ( )N is the number of nodes with degree . This system effectively keeps track of degree and heterogeneity in it, but mixing between nodes of different degrees happens at random but proportionally to degree [23, 14] , with clustering (the tendency of nodes to form connected triples) playing no role.
The degree-based mean field model can be derived exactly, under the assumption that individuals with degree k reselect their partners very rapidly, so at every moment, the status of a node is independent of the status of its current partners. In reality, many edges are long-lasting, and so correlations build up: a randomly selected infected node is more likely to connect to another infected node than a randomly selected susceptible node. More complex models are needed to correct this.
In the heterogeneous pairwise model, we also count pairs: for example [A k B ] is the expected number of links connecting a node of degree k in state A to a node of degree in state B [13, 14] ; likewise for triples of the form [A k B C m ]. The closure is done at the level of pairs (i.e. triples are approximated by singles and pairs), and hence an approximation for the triples are needed. These are given by
where ϕ is the global clustering coefficient in the network. For the un-clustered case we simply set ϕ = 0. The resulting ODEs are,
where triples are closed using equation (3) . The system can be significantly simplified for the ϕ = 0 case. When ϕ > 0, the closures become more complicated and present further challenges when implemented numerically (see notes in Appendix 6.1).
The number of equations in the heterogeneous pairwise model grows very large if the network has degrees of many different types (e.g., because there is an equation for[S k I ] for every k, pair). Generalizing this to more complex structures can become unwieldy. Edge-based compartmental models provide an alternative and are discussed next.
It is interesting to consider models which explicitly distinguish between links that happen within the households and those that happen elsewhere, as lockdowns act mostly on inter-household contacts.
To consider household structures, we take advantage of the edge-based compartmental modelling (EBCM) framework [16, 18] , adapting the model in [29] to build a model that (i) has a more realistic contact structure with households, and (ii) can distinguish between within household and community transmission, see figure 2 . This model keeps the number of equations tractable.
We assume that individuals are divided into households of size 4. Within households, there is complete mixing. In addition, each individual has a number of contacts outside the household, which allow for community transmission. The equations for this model are given in Appendix 6.2. Figure 2 : Caricature of the network model with households of size four and community stubs due to be connected up following the configuration model.
We use a version of the SEIRD compartmental model by Alvarez et al. [2] adapted to remove any built-in control measures (originally modelled as a Hill repression function modulating the number of daily contacts in response to control measures) and to include age-structured interactions in the population. The model is as follows:
where β, γ E and γ I are age-independent parameters denoting infectivity, rate at which exposed individuals become infected (inverse of incubation period) and rate at which infected individuals recover or die (inverse of disease duration) respectively (note that γ I in this model corresponds to γ in the above models). The m i are age-dependent mortality probabilities and control the fraction of those infected individuals who die. Susceptible individuals become exposed proportionally to a force of infection defined as the product of the contagion matrix with the prevalence by age. The contagion matrix is simply the product of the intrinsic infectivity of the epidemic and the daily contacts of individuals in age group i with individuals from age group j, C ij . Finally, n is the number of age groups considered and N i is the count of individuals in age group i.
Ling et al. [15] reported that the median time from symptoms onset to first negative RT-PCR in oropharyngeal swabs of convalescent patients was around 10 days, and further evidence [30] suggested that pre-symptomatic infection could happen 1−3 days before the first symptoms appear. Accordingly, we set γ = 1/14 (i.e. an average of two weeks before recovery). Before setting τ , we summarise the expression for R 0 for the various models that we consider. For the heterogeneous degree-based mean-field model [14] we have
For the heterogeneous pairwise model, we use
For the edge-based model, we set the in-household infection parameter β h to be 3 − 5 times bigger than the inter-households infection parameter β c .
The basic reproduction number of the edge-based and the age structured models are given by the leading eigenvalues of the following two next-generation matrices [8] . First, for the edge-based compartmental model, based on [24] , we have
withD and D the distribution of the excess degree (i.e. the distribution of the left-over edges attached to a node reached by following one random edge, see [22] ) and degree distribution of the network, respectively. µ 0 , µ 1 , µ 2 and µ 3 are the expected number of infected in generation 0, 1, 2 and 3 in a household of size 4. These values can be found in [4] on page 222/223, with
Finally, the R 0 for the age-structured compartmental model is given by the largest eigenvalue of:
where β is the intrinsic infectivity, γ I is the rate at which infected individuals either recover or die, N i is the size of the population in age group i, and n = 18 is the number of age groups in the model. C is the age-mixing matrix and the normalisation factor N i N j comes from the fact that at t 0 , there are only susceptible individuals in the population of each age group and therefore S i = N i in the partial derivative with respect to I of the r.h.s of the second equation in system (5).
Most of our scenarios are concerned with determining the impact of model and demographic heterogeneities on the DIHI levels. In well-mixed homogeneous populations, each individual contributes equally to spreading, and therefore DIHI is a well-defined quantity, independent of whom has been infected during the first wave. In models with degree heterogeneities, nodes with higher degree contribute to the spreading of the disease much more than nodes with fewer links. This means that depending on whom has been infected, we can observe different levels of DIHI. This effect, however, does not show in our results since the DIHI is based on the infections accumulated during the first wave as determined by our equations.
It is well known that in the simple compartmental model herd immunity at time t = 0 is achieved as long as at least 1 − 1/R 0 of the susceptible individuals are removed or immunised. In line with [3] , our general setup is that an initial epidemic spreading freely for a short time is intervened upon by implementing a lockdown period of fixed duration. Afterwards, all parameters immediately return to their pre-lockdown values. In the most basic case this is done by keeping the network or the mixing matrix the same and multiplying the transmission rate by a constant 0 < α < 1. Crucially, however, we also explore the implications of how the lockdown is modelled, that is, we investigate the difference between reducing the transmission rate whilst keeping the network the same and changing the contact network, the latter being more in tune with what happens in reality. In the more realistic age-structured model, we compare a reduction of all entries in the mixing matrix with a number of scenarios involving school closure and work distancing.
For the edge-based compartmental model we focus on final epidemic size but still under the assumption of a lockdown period. This is because we want to compare how the two different strategies affect the eventual outcome of the epidemic, rather than how the optimal α varies between the two strategies, a comparison that would be difficult to interpret.
There is an extremely strong relation between the speed or rate of spread of the uncontrolled epidemic and the timing and length of the lockdown. In fact this can be visualised in terms of the 'flattening of the curve' argument. A reduction of the transmission rate during lockdown leads to a flattening of the epidemic curve with a reduced peak and an extended duration which ideally should fit within the control period. This means that if the epidemic grows quickly and the lockdown period is short two outcomes are possible. First, a fast growing epidemic with a short lockdown period needs to be met with a significant reduction in the transmission rate, i.e., small values of α. This will lead to a reduced epidemic which does not have enough time to unfold and the lifting of lockdown is followed by a full blown epidemic. Second, if the reduction is not strong enough (i.e., larger values of α), then a significant epidemic will occur during the lockdown itself with no further peak after lifting control [7] , see also figure 3.
The main focus here is to investigate the impact of degree-heterogeneity and clustering on herd immunity induced by a first wave of the epidemic. In networks with heterogeneous degrees, one can preferentially target high-risk individuals and this will reduce the number needed to achieve DIHI. Such heterogeneity can be exploited in many ways: for example, targeted immunisations [1] and acquaintance immunisation [6] . Equally, the epidemic itself typically finds the high-risk groups first and thus 'removes' important individuals or risk groups. In line with [3, 10] , we exploit this fact and consider different levels of degree-heterogeneity and different mean-field models to explore what happens in the wake of the lockdown period when some level of spreading is possible.
In each scenario considered, we seeded nodes with degree k = 10 with [I k ](0) = 5 infected individuals, the rest of the population being fully susceptible. We let the epidemic run until the cumulative number of infected people reached 0.5% of the population. Then, a one-shot intervention lasting exactly T = 130 days kicked in. During lockdown, control measures made τ →τ 0 = ατ (by acting on τ ). Afterwards, lockdown was lifted and τ immediately returned to its pre-lockdown value. For the edge-based-compartmental and age-structured models, lockdown also involved preferential interventions on community or household links and modulation of the mixing matrix, respectively. Figure 3 shows results from the degree-based mean field model (left column) and heterogeneous pairwise model without clustering (right column) for networks with increasing levels of degree heterogeneity (from top to bottom). In each case, we find the optimal α (a simple down scaling of the transmission rate without change to the network) and report the number of infections required to achieve DIHI (i.e., total of infected and recovered nodes at the end of lockdown such that the the epidemic after lockdown is subcritical). Several observations can be made. First, for both models, aggressive control (value of α) leads to a second wave. Equally, if the control is too weak (high value of α) the epidemic will still run its course during the first wave with some reduction in the final size. Hence, there is an optimal value of α for which the final epidemic size is smallest and the epidemic post-lockdown is subcritical.
Both models clearly show that DIHI decreases with variance. Despite displaying the same overall trends, the two models are quantitatively different. In a like-for-like comparison, the degree-based heterogeneous mean-field model leads to larger overall epidemics. This is to be expected since this model does not keep track of the links explicitly and thus over estimates what would happen in a true simulation on an explicit network [14] . The pairwise model, however, accounts for links and correlations and leads to epidemics that are typically less potent. Furthermore, the heterogeneous pairwise model leads to smaller values of DIHI showing that the accuracy with which the network structure is accounted for matters. This demonstrates that model choice is important as the precise levels of DIHI matter in a real-world scenario.
The effect of clustering is illustrated by figure 4 . Typically, clustering lengthens the duration of the epidemic and lowers the peak when compared to the unclustered case (see also [29] ). The final epidemic size is also smaller. This means that there is a small amount of leeway for implementing control and that the control effort can be smaller compared to the unclustered case. It is worth noting that the final epidemic size is also smallest at the optimal α value (see also [3] ).
Finally, opting for the more accurate heterogeneous pairwise model, the level of DIHI is plotted for increasing values of variance and for different clustering levels, see figure 4 . It is clear that higher variance can drive DIHI levels to as low as 30%. Clustering in the heterogeneous pairwise model leads to even smaller values of the DIHI, although increasing variance in the degree negates the effect of clustering and the DIHI levels are very similar to those observed with ϕ = 0.25 and ϕ = 0.5. This highlights the non-trivial interactions between network properties where clustering has biggest impact in sparse networks and where high levels of degree heterogeneity can negate the effect of clustering. 
Although the previous models do account for some important contact network features, they are not ideal to capture structure such as households. Being able to capture households explicitly and having the flexibility to differentiate between household and community transmission is important because many of the interventions available to us (e.g., closing schools and workplaces) affect community transmission differently from household transmission. Thus, a distinctive feature of most lockdown measures is a change in network structure, rather than a global reduction in transmission rate.
Although our model is not an exact reflection of true population structure, it allows us to investigate whether an intervention that disproportionately affects between-household transmission can be appropriately captured by a model that treats the intervention as reducing all transmission rates.
For a given start time and intervention strength, if the duration of the intervention is long enough, the number of infections becomes very small and the eventual rebound has the same shape, regardless of duration. The optimal intervention strength leaves the population sitting at the DIHI threshold when the epidemic dies out. When the intervention is lifted, no rebound occurs. Thus the outcomes of the optimal intervention are the same if the duration is long enough. In figure 5 and at the top right of figure 6 we notice that if the intervention is long enough, the optimal final size becomes independent of duration. To see this for earlier start times requires longer durations. Our argument to explain why this happens is as follows: suppose that for a given α, we find the smallest time such that DIHI is achieved. Prolonging the intervention for longer than said time will result in the number of actively infectious people reducing even further by end of lockdown but not change the downward trend of the epidemic, which will eventually die out.
Comparing the top right of both panels in Figure 6 show that the optimal community intervention allows more infections than the optimal global intervention.
The intuition behind this is based on the observation that epidemics typically exploit 'heterogeneities' in the population. For networks,, this means that high-degree nodes typically become infected early on in the epidemic. This is the main reason why a first epidemic wave in a network with very heterogeneous degree infects the highly connected nodes. In a scale-free like network, the number of such nodes is small (e.g., 20% of the nodes responsible for 80% of infection). If the epidemic progresses with strong interventions in place, it cannot spread far beyond these high-degree nodes. Once it dies out, the residual network is highly fragmented and made up of much lower degree nodes.
In this household model, community links drive degree heterogeneity. Household links alone lead to a regular network. Hence, when we effectively cut most community links, heterogeneity in degree is significantly reduced for the duration of control. This means that many high-degree nodes that would normally be infected during the first wave will now not get infected. The infection is unable to target the highest-degree nodes. When control is lifted, the high-degree nodes reactivate their links, allowing the epidemic to rebound.
However, when controlling both link-types equally, degree heterogeneity is preserved and the infection again preferentially targets the high-degree nodes. A weaker control targeted to all edges rather than just community edges may allow an initial wave of similar size, but it will preferentially target the high-degree nodes. So this type of control acts as a very effective way of finding highly-connected nodes. This then means that at the end of lockdown, it is more likely that the most 'dangerous' nodes be removed, and with that, a smaller chance of a second wave.
Consider two epidemics with the same intervention start time, providing the optimal strategy for either the global or the community case. We would expect that the average community degree of those who have been infected in the global case would be higher than in the community case. Thus, on average, the individuals immunised by infection in the community case will be less important to disease transmission and thus more of them must be immunised to achieve the DIHI threshold. This is important because households may be able to sustain the epidemic for extended amounts of time and therefore change the outcome in DIHI levels.
In this section, we further explore the notion that modulating the effective R 0 of the epidemic through modifying the structure of the mixing matrix (in this case, the age-structure mixing matrix) can affect the system differently from achieving it through simply scaling each element of the mixing matrix. To do so, we began by considering three scenarios described by [26] , namely, school closure, school closure and social distancing, and work distancing. In a baseline, no-intervention case, the matrix of daily contacts C was set to be the sum of 4 components: school contacts, work contacts, home contacts and other contacts. In what follows, we use the corresponding matrices of agebanded daily contacts in the UK produced by the POLYMOD study [26] . Briefly, school closure is realised by zeroing the school component of the mixing matrix); school closure and social distancing involved zeroing school contacts as well as reducing by half the number contacts at other locations between school-going individuals (first four age groups); work distancing is implemented by halving the contacts made at the workplace. For each of these interventions, we compared the behaviour of the system when re-scaling the matrix of total contacts so that its R 0 during the intervention was the same as the R 0 of the modified mixing matrix during the intervention (namely, 2.114, 2.106 and 2.179 for scenarios 1 to 3 respectively).
Age-banded population counts (18 5-year bands) were taken from the Office for National Statistics (licensed under the Open Government Licence), pooling all age groups above 85. Mortality rates were taken from the modelling of [28] , assuming that the rates in ten year age bands are the same as across two five year age bands. These rates were calculated based on cases in China in the initial outbreak and from the closed population on the Diamond Princess cruise ship. It is possible that these rates may prove to be overestimates compared to populations that do not experience overwhelming levels of hospitalisation, but they are not likely to be impacted by the interventions considered here.
To determine the evolution of the epidemic in this baseline case, we first scaled the contact matrices so that the system's R 0 was 2.5 (to maintain consistency with previous sections). In all cases, we used the same start date (T = 260) and duration (150 days) for the intervention. These parameters were arbitrarily chosen among the sets of possible parameters resulting in a sub-critical epidemic post-intervention. Figure 7 confirms that all interventions (coloured lines) result in a reduction in the number of infected individuals. However, whether this intervention is realised through modulating the contact matrix (solid lines) or through uniform scaling (dashed lines) results in substantially different outcomes (≈ 16% differences for scenarios 1 and 2, ≈ 3% in scenario 3), even though the scaling factors used in the control interventions (dashed lines) are very similar from one scenario to the other. This is clear evidence that the structure of the contacts modulates the effects of the intervention. Figure 8 provides a visual intuition as to why the work distancing scenario is closest to simply scaling the matrix. As pointed out by [21] , assortative mixing dominates in 3 of the components (home, school, other). Thus, school closure primarily affects diagonal elements of the mixing matrix (and primarily for the first 4 age groups). Intervention 2 does involve halving (some of) the contacts in the other component and some of those terms diffuse away from the diagonal, however, these contribute little to the overall mixing matrix. In contrast, the work component is the only component to feature what Mossong et al. [21] describe as a wide contact plateau. Because this plateau accounts for more than half of the total number of contacts within the corresponding 8 age-bands, intervention 3 (social distancing) is most akin to scaling the entire matrix.
To further clarify how the structure of the contacts modulates the effects of the intervention, we carried out simulations in which two confounding factors were removed, namely heterogeneity in the number of individuals in the different age-bands and in the frequency of contacts by age. Whereas the former plays a key role in the calculation of the effective R 0 (see Section 2.6), the latter weights the impact of the intervention. For example, zeroing school contacts which only account for 12% of the total number of contacts will be negligible compared to zeroing other contacts that account for 40%. Therefore, in what follows, all age groups were set to have the same number of individuals (1/18-th of the total population) and all contact components were scaled to have the same sum of elements (arbitrarily, the sum of the original number of other contacts). We then systematically analysed the effect of three different interventions in which one component (school, work, or other was systematically scaled down by a factor taking values between 1.0 (no intervention) and 0.0 in steps of 0.1. Figure 9 shows the effect of the most severe form of intervention (the zeroing of the relevant component) and clearly demonstrates that, once confounding factors are removed, changes in the mixing matrix lead to different outcomes in terms of whether such an intervention is more or less effective than simply scaling the matrix to achieve the same effective R 0 . Here, zeroing the school component is less effective than scaling the overall matrix (despite an overall scaling factor of 0.61). Instead, zeroing either the work or other component is more effective than scaling the overall matrix (despite larger scaling factors, 0.87 and 0.88, respectively). It should be noted that whilst zeroing the work and other components leads to similar results in terms of total prevalence, there are differences in prevalence by age which could have significant implications when age-structured mortality rates are considered.
In this paper we explored a range of mean-field models previously used to approximate exact epidemics on networks and providing some analytical traction regarding how network properties impact epidemic invasion, final size and the efficacy of control measures. In order of increasing complexity, these models are: the degree-based heterogeneous mean-field, heterogeneous pairwise (without and with clustering), and an edge-based compartmental model which explicitly includes household structure and can distinguish between household and community transmission. While these cannot be used as they are to inform policy they still provide important insights into model selection and key features that need to be captured, or can be exploited, to identify the best possible control measures. In addition, we also tested our findings against a more realistic age-structured model with real mixing matrices.
We have shown that increased degree heterogeneity (i.e., higher variance in the degree distribution) leads to DIHI levels that are much smaller than the basic compartmental model, 1 − 1/R 0 . This is line with the findings of [3, 10] . Moreover, we quantified the extent to which the DIHI induced by the first wave depends on the variance in the degree distribution. We have shown that herd immunity in clustered networks is even lower because epidemics on clustered networks last longer and have lower peaks, allowing more flexibility regarding the start and intensity of control.
Perhaps, the most important question that we addressed regards how lockdown/control is implemented in different models. Many models assume that during lockdown the contact network or mixing matrix is not changing but rather the transmission rate is scaled [3, 7, 10, 20] . We do not believe that this is appropriate because during lockdown the underlying contact structure changes. Our results with the edge-based compartmental and age-structured models have shown that these two approaches differ in outcome. Perhaps, the assumption of a non-changing contact structure during lockdown is more likely to be made in mean-field models. In models at higher than mean-field resolution (e.g. agent-based) it is much easier to explicitly modify the contact network.
We have therefore revealed a possible risk with using a model that ignores household structure to infer the level of infection required to reach the DIHI threshold. The favourable change in the DIHI threshold compared to what we predict from a homogeneous model is a consequence of heterogeneities. Where the intervention makes the population more homogeneous, the disease will no longer act like an effective intervention. In an extreme case, it will require more infections to achieve herd immunity than a random vaccination would need [9] . As the infection spreads along edges, we would see that at the end of the first wave, susceptible people are disproportionately in contact with other susceptible people and recovered people are disproportionately in contact with recovered people. So the residual network of susceptible nodes has more contacts than would occur if the same fraction were vaccinated randomly.
Another important observation resulting from our work is that it is extremely difficult to make general statements by extrapolating from findings based on simple models. Most models in fact ignore meso-scale structures (e.g., degree heterogeneity does well for local or micro structure whilst mixing matrices do well for macro-scale mixing) and their absence may exacerbate the impact of an intervention (either positively or negatively) leading to erroneous conclusions. In the present paper we saw that when intervention could not act on the global network of contacts, DIHI levels varied substantially, although heterogeneities still played a major role in reducing them.
Finally, it is worth noting that in all models we considered, when lockdown ended (and if DIHI was not reached), epidemics went on to grow exponentially. However, in many real-life scenarios, a prolonged phase of sub-critical spreading (i.e. slow decay) has been observed before the exponential growing phase returned. There is a number of reasons why this behaviour is not observed in our models: (i) Deterministic models fail to capture fluctuations that dominate when the number of infected people is small, which might have a major impact on resurgence; (ii) during lockdown, the contact structure changes drastically and abruptly, but when lockdown is lifted, there is a delay/inertia in going back to pre-lockdown status, such lag having important implications on resurgence; (iii) after lockdown, some social distancing measures remain along with social awareness reducing exposure to the disease. Understanding and modelling these effects is an interesting challenge to address in future work.
In complementing our study of network-based mean-field models, we used a model with no explicit contact structure, where instead contact structure was alluded to via age-related mixing patterns. In the present situation with Covid-19, such models are appealing because they explain some of the structure of the population and provide a rationale for establishing a lockdown. Moreover, with higher mortality rates among older people, age-structured models are of interest in their own right. There is some tension between which model is most apt for describing population experience of the infection, the control measures effected and the outcomes for the population. The ideal solution may lie between some of the options presented here.
However, age-structured models cannot say much that is explicit about the structure of contacts within the population. While we have discussed ways in which lockdowns can be implemented in such models, the formulation is arguably less intuitive than in the network case. How to bridge the gap is an interesting question. An age-related network might be one labelled with age classes, with analysis focused on understanding the positions in the network occupied by individuals of a given age. This could be coupled to household models of network formation, one with variable sizes and smaller households more likely to feature older individuals. Once the general structure of such networks is known, adapting the models here would be simple enough. These considerations are beyond the presentation here but would be a fruitful avenue for future discussions, particular if additional lockdowns are required.
We can account for clustering in the pairwise model, by introducing a clustering factor ϕ. The equations are similar to (4), but when closures are implemented we have to consider both open and closed triples.
In the case of open triples, the equations for the closures are the same as in (4), scaled by (1 − ϕ) to account for clustering, i.e.
where the superscript o indicates open triples. For closed triples, we first study [A S k I m ] c . This quantity represents triples in which the first node is in state A, the middle node is in state S with degree k, and the third node is in state I with degree m. The triple is closed, therefore the third node is connected to the first one. This introduces correlations when we write the triple in terms of pairs for the closure, which we indicate with C S I M .
[
The correlation can be written in terms of the ratio between realized pairs [A k I m ] and possible pairs [A k I m ] in a well-mixed population:
Hence, the closure is
In a similar manner, we can write the expression for [I m S k A ]. The resulting system is thereforė 
We consider a network with 4N nodes partitioned into households of size 4. Apart from the within household community, each node has a number of links to nodes outside the household, according to the degree distribution P n,p (k). The within household per-contact-transmission is denoted by β h while the community transmission by β c . The resulting system is given below,
ϕ RRR (t) = γϕ IRR (t),
ϕ SII (t) = −(A(t) + 2β h )ϕ SII (t) − 2β h ϕ SII (t) + 2(A(t) + β h )ϕ SSI (t) − 2γϕ SII (t), (24) Θ(t) = −(β h ϕ SSI (t) + 2β h ϕ SII (t) + 3β h ϕ III (t) + β h ϕ SIR (t) + 2β h ϕ IIR (t) + +β h ϕ IRR (t)),
A(t) = ψ (θ(t)) ψ(θ(t)) β c ϕ I (t),
S(t) = (1 − )Θ(t)ψ(θ(t)),
with the following initial conditions:
with all other variables set to zero at time t = 0. 