The importance of saturating density dependence for population-level predictions of SARS-CoV-2 resurgence compared with density-independent or linearly density-dependent models, England, 23 March to 31 July 2020

Background Population-level mathematical models of outbreaks typically assume that disease transmission is not impacted by population density (‘frequency-dependent’) or that it increases linearly with density (‘density-dependent’). Aim We sought evidence for the role of population density in SARS-CoV-2 transmission. Methods Using COVID-19-associated mortality data from England, we fitted multiple functional forms linking density with transmission. We projected forwards beyond lockdown to ascertain the consequences of different functional forms on infection resurgence. Results COVID-19-associated mortality data from England show evidence of increasing with population density until a saturating level, after adjusting for local age distribution, deprivation, proportion of ethnic minority population and proportion of key workers among the working population. Projections from a mathematical model that accounts for this observation deviate markedly from the current status quo for SARS-CoV-2 models which either assume linearity between density and transmission (30% of models) or no relationship at all (70%). Respectively, these classical model structures over- and underestimate the delay in infection resurgence following the release of lockdown. Conclusion Identifying saturation points for given populations and including transmission terms that account for this feature will improve model accuracy and utility for the current and future pandemics.


Introduction
Like many pathogens that cause respiratory diseases [1][2][3], severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) appears to be transmitted more effectively in densely populated areas [4][5][6]. The increased disease rates reported among high-density populations [4,5,7,8] may, however, be an artefact of confounders, such as a higher proportion of individuals of lower socioeconomic status or from minority ethnic groups in urban areas [9]. Using coronavirus disease (COVID-19)associated mortality data from the Office for National Statistics, we aimed to assess the evidence for density dependence.
Standard transmission models that either do or do not account for this density dependence have been used interchangeably because their projections are generally equivalent when population density remains unperturbed or is homogeneous, e.g. at a national level. While the ca 1% infection fatality rate for COVID-19 [10] is insufficient to destabilise populations, the reaction of most countries' governments to curtail disease spread through lockdown and physical distancing has had unprecedented impacts on the density of mobile human populations. For example, the United Kingdom's lockdown, which came into effect on 23 March 2020, effectively reduced the freely moving population from 66.5 million to 10.6 million (key workers) [11]. This same intervention was employed by numerous countries, similarly impacting their mobile populations [12]. We evaluate the extent to which models built to inform the epidemiology of COVID-19 use an underlying structure that can accommodate the drastic changes and variation in densities experienced by most global populations.
As lockdowns were gradually released over the latter part of 2020, global populations were expected to re-equilibrate to a 'new normal' whereby densities of mobile people were increased but in which contact patterns were expected to remain reduced through physical distancing interventions [13]. Using a suite of mathematical models, we illustrate the impact that the different, routinely ignored, assumptions underlying transmission and density may have in projecting infection dynamics and measuring intervention effectiveness.

Data
Reported COVID-19-related deaths between 1 March and 31 July 2020 were obtained in anonymised linelist form from Public Health England and were filtered to include all deaths which occurred within 28 days of positive COVID-19 test (n = 36,311). We aggregated individual records to lower-tier local authority (LTLA), and nationally by 10-year age bands in order to calculate age-standardised expected counts.
Local authority shapefiles and single-age population estimates were obtained from the Office for National Statistics [14]. Four sub-regions of Buckinghamshire (Aylesbury Vale, Chiltern, South Bucks, Wycombe) were aggregated in order to match most recent population estimates. The City of London was aggregated with Westminster because of its very small resident population, and the Isles of Scilly were excluded since no COVID-19-related deaths had been reported there during the study period. Index of multiple deprivation (IMD) [15], percentage of minority ethnic population [16] and percentage of key workers among the working population [17] are characteristics of the LTLA population potentially associated with both COVID-19 mortality and population density, therefore we included them as covariates in all models. Percentage of key workers was missing for Westminster and Cornwall; these were imputed by the median value across all neighbouring LTLA.

Statistical analysis
Negative binomial regression models were fitted to the number of deaths (n) per LTLA, adjusting initially for the three covariates (IMD, % minority population, % key workers) and subsequently adding a fourth covariate, namely the lag in weeks behind the first death nationally. We adjusted for age distribution within the LTLA via inclusion of age-adjusted expected deaths (E) as an offset; these were calculated according to national age-specific rates (deaths per 100,000 per age band) applied to local population estimates in 10-year age bands. We accounted for population density in one of four functional forms: (i) constant/independent of population density, (ii) linear, (iii) log-linear and (iv) saturating.
For observed number of deaths (D), age-adjusted expected deaths (E) and defining x i IMD , x i mino , x i KW , x i lag and x i dens , respectively, as the deprivation score, % minority population, % key workers, lag in weeks behind first death nationally and population density of LTLA this yields the following model specification: where μ i is the rate relative to the expected, ϕ the size parameter (1/overdispersion) of the negative binomial distribution and β j the regression coefficients. The functional form of population density is defined as where Model compartments are: 'S'usceptible, 'E'xposed, 'I'nfectious, 'P're-critical infectious, 'C'ritically ill, 'D'ead and 'R'ecovered. The three functional forms are shown for the force of infection, λ. Each is a product of the infectious proportion of the population and the transmission coefficient, β. The transmission coefficients are differentiated between the three functional forms using primes. The transmission coefficient is itself the product of the transmission probability per contact with an infectious individual and the contact rate. Under a linearly density-dependent assumption, the contact rate increases linearly with the total population, N. This numerator N cancels out the N denominator which is why β` is multiplied by the number of infectious individuals instead of the proportion. The saturating density-dependent formulation assumes the force of infection is a product of the transmission coefficient and a function that increases contact rate nonlinearly with population density as informed by parameter θ derived from analysing England's regional mortality data.
Models were fitted using the rstanarm package [18] with default weakly-informative priors. We compared the four model variants on leave-one-out information criterion (LOOIC), calculated via approximate leave-one-out cross-validation as implemented in the loo package [19]. Interpretation was the same as that of the Akaike information criterion in that smaller values reflected better fit. The value of ϴ for the saturating function was determined by manual optimisation of the generalised linear model with respect to LOOIC on a hold-out set of 40% of LTLA, over a range from 0.001 to 1.
For the saturating model, the impact of an 84% reduction in effective population density as a result of lockdown on predicted mortality rates among the freely moving population was calculated as a percentage change between mean model-predicted deaths under the original and reduced densities.

Mathematical model
We use a discrete-time, deterministic compartmental model ( Figure 1) with daily timesteps to simulate SARS-CoV-2 transmission. From the first day of lockdown (23 March 2020), we assumed 84% of the population to enter isolation in which frequency-dependent transmission occurs. We made this assumption for the lockdown sub-population because an individual's likelihood of contracting infection while in their home is limited by their household size (i.e. not impacted by the density of individuals under isolation in different households). Each model was fitted independently to England's COVID-19-associated mortality data (up until 1 August 2020). We compared frequency-dependent and both linearly and saturating density-dependent transmission for when lockdown is released. We also explored the impact of varying rates of connectivity between locked-down and free-moving individuals because those under lockdown were still afforded some freedom of movement, and because key workers potentially cohabit with those under lockdown (Supplementary Figure S1 shows limited impact of doubling the rate of this connectivity). We compared frequency-dependent and linearly density-dependent transmission (the limiting cases for the saturating density-dependent model [20]) among the remaining free-movers for a range of lockdown release schedules (over a period of between 1 and 12 months). Contact rates are reduced through two distinct mechanisms under the density-dependent models: whereas reduced contact through physical distancing and behavioural changes among the freely moving population (e.g. the 2 m rule) was included in all models, only the densitydependent versions assumed reduced opportunities for mobile people coming into contact with other mobile people because of their substantially depleted numbers. Full model specification and sources for its parameterisation can be found in the Supplement, and the Python (v3.8) code is freely available from github (https://github.com/lwyakob/COVIDsaturates).

Ethical statement
No ethical clearance was needed for this publication because all data which were freely available from the Office for National Statistics were anonymised.

Results
Evidence for saturating density dependence in COVID-19-associated deaths Owing to the heightened risk of earlier outbreak seeding for higher density areas, we repeated the analysis, additionally adjusting for the lag of the local epidemic behind the national. The saturating model was retained as the best fit (Table) and suggested a similar increase in mortality rate of fourfold (90% CrI: 2.15-7.08).
Under the saturating density-dependent model, the impact of lockdown on reducing transmission among mobile individuals, and consequently on deaths, is heterogeneous, having greatest benefit in regions with low population density (> 30% reduction in projected deaths for example in Devon, Herefordshire and the Derbyshire Dales) but reduced benefit in high-density regions (ca 5-7% reduction for the London boroughs of Tower Hamlets, Hackney, Islington and Camden) ( Figure 2B). These results were retained when accounting for the lag of the local epidemic behind the national (Supplementary Figure S3).

Projecting SARS-CoV-2 resurgence after lockdown is released
A full-text review of 100 epidemiological models of SARS-CoV-2 published until 19 June 2020 showed that 70% explicitly assumed that contact rate between people (and hence transmission) is unaffected by population density (see Supplement for details of the models). Of the remaining 30% of models, all assumed a linear relationship between population density and transmission.
We used a metapopulation model to simulate the infection dynamics among freely moving as well as lockeddown individuals, incorporating transmission terms that can accommodate density-independent (referred to as 'frequency-dependent') as well as linearly and saturating density-dependent assumptions. While all functional forms performed equivalently in fitting mortality data leading up to lockdown, dynamics under alternate assumptions may diverge markedly during and following the phase when lockdown is released ( Figure 3). We note that we ignored any adaptive public health responses (i.e. additional interventions) curbing the second wave -this comparison was intended to illustrate the consequences to projected dynamics of alternative assumptions underlying density and transmission.
Although final epidemic size and total deaths were equivalent for the alternative classical assumptions (fd and linearly dd), transmission was delayed by almost a year under a density-vs frequency-dependent model ( Figure 3). This delay occurs because only under the density-dependent assumption, the force of infection is reduced while any part of the population remains locked down. At the very high densities of London populations, locking down 84% of people under our saturating density-dependent model had an impact most similar to a frequency-dependent assumption. This means that, if the density of England's entire population was equivalent to the density found in London, infection dynamics and deaths resulting from a saturating density-dependent model most closely match the frequency-dependent projections. However, London has a population density that is an order of magnitude higher than the next most populated region in England, and projected infection dynamics diverged more considerably under scenarios reflecting densities experienced outside of the capital. The force of infection and the timing of peak prevalence for the saturating density-dependent model is constrained between the frequency-and linearly density-dependent versions [20], with lower densities tending towards the latter.
Assuming a maximum national capacity of 5,000 intensive care unit (ICU) beds, we assessed the difference between these temporal limits in the projected duration between the release of lockdown and a second wave of infection exceeding ICU capacity (Figure 4). Threshold levels of physical distancing to interrupt transmission (i.e. maintain an effective reproduction number < 1) are similar for both classical models. However, where interventions failed to achieve this threshold, density-dependent transmission resulted in a delay of more than a year before the ICU capacity was exceeded. This is contingent on the timeframe across which lockdown is released, whereby more gradual releases extend delays.

Discussion
Projections of COVID-19 infection dynamics following the release of a huge proportion of the population from lockdown comprise an urgent and critical component of public health decision-making [22]. The classical forms of modelling infectious diseases among populations have been used interchangeably by different research groups because, under most plausible circumstances, they exhibit equivalent dynamics. In March 2020, England locked down more than four-fifths of its population. For most, this fundamentally altered the rate at which people made contact. Under the circumstance of millions of people easing out of lockdown, substantial differences between projections from a frequency-and density-dependent transmission assumption emerge. Most notably, density dependence results in delayed infection resurgence and, contingent on the timeframe across which lockdown is released and the effectiveness of physical distancing, this delay can extend to more than a year.
The delay is a function of a fundamental aspect of density-dependent transmission: lower host densities reduce the force of infection, and there is a threshold host density below which an infection cannot spread. Despite its origins in human infectious disease modelling [23], the existence of this threshold has historically had limited epidemiological application. The phenomenon is discussed more widely in wildlife disease ecology [24] where it underlies key disease control decisions such as culling [25]. Current expectation is that lockdowns, either full or in a more moderate or localised form, will be reimplemented when cases start increasing again. Density effects and thresholds are particularly pertinent in the current COVID-19 pandemic during which extreme fluctuations in mobile human density are likely to continue.
Analysing COVID-19-associated deaths across different regions in England, and accounting for known major confounders [9], the nonlinear increase in deaths with population density was adequately captured by neither classical form of modelling transmission. Using a function that captures the saturating increase in deaths with population density resulted in an expedited resurgence compared with a linearly density-dependent model and a delayed resurgence compared with the popularly used frequency-dependent model.
Less populated areas were shown to have fewer deaths per capita (as per England's mortality data) and slower resurgences following the release of lockdown. This provides more achievable targets and considerably more lead time for health services to prepare than would otherwise be anticipated. It also highlights a hazard. During and after releases from lockdown, in order to fit a prolonged lag in cases, transmission rates derived from most current (frequency-dependent) models will underestimate the effective reproduction number. This could exaggerate the perceived effectiveness of ongoing interventions, such as physical distancing or face masks, with potentially serious consequences.
Our study is limited by the fact that we do not have comprehensive data on how contact rates were affected before and over the lockdown period for

Figure 3
Population density and SARS-CoV-2-associated mortality and infection dynamics following the release of lockdown, by transmission assumptions, England, March-July 2020 individuals inhabiting regions of differing population density. We also do not know where people were infected, only where they were when they died. Instead we had to resort to mortality rates and locations as a proxy. It is possible, for example, that contact rates are not affected by dramatic shifts in population density regardless of baseline levels (i.e. the average England resident came into contact with as many individuals during lockdown as before lockdown, satisfying a frequency-dependent assumption), and that the increased per capita fatality seen in more densely populated regions has an alternative, thus far unidentified, explanation. Mobile phone applications developed to notify participants of urgent health information have already gained millions of users in the current coronavirus context [26]. Piggy-backing on these efforts could help substantiate the evidence for the contact-density relationships we have identified.
Owing to the highly complex interactions between population characteristics, behaviours and mortality risk, the association discovered between saturating density and mortality rates may remain confounded by factors not considered here. Moreover, the criterion used for model comparison depends on an independence assumption which may not hold between neighbouring LTLA. Work is ongoing to characterise the patterns of spatial correlation in mortality at the LTLA level.

Conclusion
Infectious diseases are emerging at an unprecedented rate and the upwards trend in global travel and urbanisation increases the likelihood of pandemics. Their success in controlling SARS-CoV-2 means that widescale lockdowns will not only continue to be enforced as this pandemic progresses, but they are likely to be more readily applied in future emergencies. It is crucial that we use the current opportunity to collect data to inform more precise forms of how contact rates are altered at varying stages of lockdown. Future work should also address whether the feature of saturating density dependence we have identified from England's data are generalisable to other countries. Incorporating realistic contact-density relationships into the transmission term of population-level mathematical models will improve precision of their projections and their utility in public health decision-making.

Figure 4
Consequences of density dependence on intensive care unit capacity, England, March-July 2020 Difference (in days, contours) between frequency-and densitydependent models in their projected time before ICU capacity is exceeded, as a function of lockdown release schedule and effectiveness of physical distancing. The dashed vertical line marks minimum physical distancing required to prevent exceeding ICU capacity in the frequency-dependent model.