--- title: "Limitations of models for guiding policy in the COVID-19 pandemic" author: - name: Paul M McKeigue affiliation: Usher corresponding: yes email: paul.mckeigue@ed.ac.uk - name: Simon N Wood affiliation: Math address: - code: Usher address: Usher Institute, College of Medicine and Veterinary Medicine, University of Edinburgh, Teviot Place, Edinburgh EH8 9AG, Scotland, UK - code: Math address: School of Mathematics, James Clerk Maxwell Building, University of Edinburgh, EH9 3FDG, Scotland, UK header-includes: \usepackage{newunicodechar} \usepackage[T1]{fontenc} \let\origquote\quote \def\quote{\origquote\itshape} \usepackage{graphicx} \usepackage{geometry} \usepackage{longtable} \usepackage{booktabs} \usepackage{float} \floatplacement{figure}{H} \usepackage{array} \usepackage{threeparttable} \usepackage{longtable} \usepackage{setspace} \usepackage{soul} % \usepackage{newtxtext,newtxmath} # \renewcommand\hl[1]{#1} #\newcommand{\pdif}[2]{\frac{\partial #1}{\partial #2}} output: rticles::plos_article: includes: # in_header: ./preamble.tex latex_engine: lualatex keep_tex: true fig_caption: yes md_extensions: -autolink_bare_uris geometry: "left=2cm,right=2cm,top=2cm,bottom=2cm" bibliography: ./covid.bib link-citations: yes csl: ./biomed-central.csl #csl: ./vancouver.csl always_allow_html: true urlcolor: blue linenumbers: false linkcolor: cyan --- ```{r functions, message=FALSE, echo=FALSE} library(knitr) library(kableExtra) library(data.table) library(ggplot2) options(knitr.kable.NA = '.') options(knitr.table.format = "latex") # adds tab: prefix to labels knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE) knitr::opts_chunk$set(fig.width=6, fig.pos = "H", out.extra = "") ## lambda must be < alpha ## -R0 (1 - S) is negative for S < 1, so this always holds mgf.gamma <- function(theta, alpha) { if(is.infinite(alpha)) { mgf <- exp(theta) } else { mgf <- (1 - theta / alpha)^-alpha } return(mgf) } finalsize <- function(R0, lambda=1) { if(lambda==1) { alpha <- Inf S.minusSinf <- function(S) (S - exp( -R0 * (1 - S)))^2 } else { alpha = 1 / (lambda - 1) S.minusSinf <- function(S) (S - (1 + R0 * (1 - S) / alpha)^-alpha)^2 } S.inf <- optim(par=0, # start from S=0 so that we find the root that is < 1 fn=S.minusSinf, lower=0, upper=1, method="Brent", control=list(trace=1))$par return(1 - S.inf) } runepidemic <- function(R0, lambda=1, rategamma=0.2, uncorr=TRUE) { numdays <- 200 dSt <- numeric(numdays) dqt <- numeric(numdays) St <- numeric(numdays) It <- numeric(numdays) Vt <- numeric(numdays) epsilon <- 0.0001 if(uncorr) { alpha <- ifelse(lambda > 1, 1 / (lambda - 1), Inf) # uncorrelated beta0 <- R0 * rategamma } else { alpha <- ifelse(lambda > 1, 2 / (lambda - 1), Inf) # connectivity beta0 <- R0 * rategamma / (1 + 1 / alpha) } H.S <- R0^(-1 / lambda) # fraction susceptible at herd immunity threshold St[1] <- 1 It[1] <- epsilon Vt[1] <- 0 propfall.beta <- 0.5 beta.t <- rep(beta0, numdays) for(i in 1:(numdays-1)) { hS <- ifelse(lambda > 1, St[i]^lambda, St[i]) if(uncorr) { dSt[i] <- - beta.t[i] * hS * It[i] dqt[i] <- - beta.t[i] * It[i] } else{ dSt[i] <- - (1 + 1 / alpha) * beta * hS * It[i] dqt[i] <- - beta.t[i] * (1 + 1 / alpha) * St[i]^(1 + 2/alpha) * It[i] } qt <- alpha * (1 - St[i]^(-1/alpha)) dVt <- rategamma * It[i] dIt <- -dSt[i] - dVt St[i+1] <- St[i] + dSt[i] It[i+1] <- It[i] + dIt Vt[i+1] <- Vt[i] + dVt } epidem.dt <- data.table(t=1:numdays, cases=-dSt, beta.t, St, dqt, qt, It, Vt, lambda=lambda, uncorr=uncorr, alpha=alpha, H.S=H.S) epidem.dt[, R := R0 * St^lambda * beta.t / beta0] epidem.dt[, H.t := t[which.min((St - H.S)^2)]] # time at which herd immunity threshold is reached epidem.dt[, lambda := factor(lambda)] return(epidem.dt) } lambda <- c(1, 4) R0 <- 2.5 rategamma <- 0.2 # generate epidemics for values of lambda epidem.dt <- NULL for(i in 1:length(lambda)) { epirun <- runepidemic(R0, lambda[i], rategamma, uncorr=TRUE) epidem.dt <- rbind(epidem.dt, epirun) } epidem.H <- unique(epidem.dt, by="H.S") ## calculate final size for sequence of values of R0 and lambda R0 <- c(1.2, 1.5, 2, 3, by=1) lambda <- seq(1, 10, by=0.1) fsize.het <- NULL for(j in 1:length(R0)) { for(i in 1:length(lambda)) { fsize.het <- rbind(fsize.het, data.table(R0=R0[j], finalsize=finalsize(R0[j], lambda=lambda[i]), lambda=lambda[i])) } } ``` # Abstract At the outset of the COVID-19 epidemic in the UK, infectious disease modellers advised the government that unless a lockdown was imposed, most of the population would be infected within a few months and critical care capacity would be overwhelmed. This paper investigates the quantitative arguments underlying these predictions, and draws lessons for future policy. The modellers assumed that within age bands all individuals were equally susceptible and equally connected, leading to predictions that more than 80% of the population would be infected in the first wave of an unmitigated epidemic. Models that relax this unrealistic assumption to allow for selective removal of the most susceptible and connected individuals predict much smaller epidemic sizes. In most European countries no more than 10% of the population was infected in the first wave, irrespective of what restrictions were imposed. The modellers assumed that about 2% of those infected would require critical care, far higher than the proportion who entered critical care in the first wave, and failed to identify the key role of nosocomial transmission in overloading health systems. Model-based forecasts that only a lockdown could suppress the epidemic relied on a survey of contact rates in 2006, with no information on the types of contact most relevant to aerosol transmission or on heterogeneity of contact rates. In future epidemics, modellers should communicate the uncertainties associated with their assumptions and data, especially when these models are used to recommend policies that have high societal costs and are hard to reverse. Recognition of the gap between models and reality also implies a need to rebalance in favour of greater reliance on rapid studies of real-world transmission, robust model criticism, and acceptance that when measurements contradict model predictions it is the model that needs to be changed. \newpage # Introduction At the outset of the COVID-19 epidemic in the UK in early 2020, advice from the Scientific Advisory Group on Emergencies (SAGE) to the UK government relied on models that predicted that unlesss coercive restrictions of social contact ("lockdown") were imposed without delay, the epidemic would infect most of the population within a few months and critical care facilities would be quickly overwhelmed [@ferguson_report_2020;@davies_effects_2020]. Together with reports from northern Italy that hospital services had been overwhelmed with cases requiring critical care, this led to an abrupt change of policy and the imposition of a lockdown starting on 24 March 2020. Subsequent management of the pandemic in the UK and other countries was based on suppressing transmission through restricting social contacts, vaccinating the entire adult population and a test and trace programme, rather than on focused protection of the vulnerable. Reversal, or even questioning, of these policies was difficult, in part because of implementation of the recommendation from SAGE's behavioural science subgroup to attempt to overcome people’s own rational assessment of risk with “hard hitting emotional messaging” [@spi-b_options_2020]. The rapid end of the first wave of the epidemic was attributed to lockdown, despite evidence by early May that infection rates had been in substantial decline before lockdown. As similar epidemic modelling approaches may be used to advise policy makers about emerging epidemics in the future, it is important to examine their limitations. The advice that only a lockdown could prevent critical care capacity from being overwhelmed in the UK, presented in reports from modellers at Imperial College (IC) [@ferguson_report_2020] and the London School of Hygiene & Tropical Medicine (LSHTM) [@davies_effects_2020] was based on four propositions: (1) in an unmitigated epidemic with basic reproduction number ${\cal R}_0=2.5$, more than 80% of the population will be infected within a few months; (2) about 2% of those infected will require critical care; (3) to limit morbidity and mortality will require suppressing the epidemic by reducing the reproduction number below 1, rather than mitigation by shielding the vulnerable; (4) only a lockdown can reduce the reproduction number below 1. In a subsequent report the IC modellers asserted that these propositions applied globally [@walker_report_2020]. This article examines the quantitative arguments underlying these propositions and the lessons that can be learned. # Methods The approach used by the IC and LSHTM modellers was based on the classic Susceptible-Exposed-Infected-Removed (SEIR) compartmental model formulated by Kermack and McKendrick [@kermack_contribution_1927] in 1927, in which the **transmission function** relating the infection rate to the proportion $S$ of individuals remaining susceptible is assumed to be linear in $S$ [@anderson_reproduction_2020;@keeling_predictions_2021]. This is equivalent to assuming that all individuals are equally susceptible and equally connected [@novozhilov_spread_2008]. Kermack and McKendrick were well aware of this, stating in the abstract of their 1927 paper: > In the present communication discussion will be limited to the case in which all members of the community are initially equally susceptible to the disease. From their model, Kermack and McKendrick derived an equation (Appendix Equation \ref{eqn:fsizehom}) for the final size of an unmitigated epidemic given the basic reproduction number ${\cal R}_0$. With ${\cal R}_0=2.5$, the herd immunity threshold is `r round(100 * (1 - 1 / 2.5))`%, and the epidemic ends when `r round(100 * finalsize(2.5, 1))`% of the population have been infected. The assumption that susceptibility and connectivity do not vary between individuals is unrealistic; relaxing this assumption to allow for heterogeneity gives rise to a **non- linear** transmission function [@novozhilov_spread_2008;@gomes_individual_2022;@neipel_powerlaw_2020;@tkachenko_timedependent_2021] relating the infection rate to the proportion $S$ remaining susceptible. This function is mathematically constrained to be non-negative, increasing, convex (second derivative non-negative), and smooth unless a discrete factor like vaccination has a large effect. Under these constraints the function can be approximated by a power law in which the infection rate is proportional to $S^\lambda$, where $\lambda$ is the **immunity coefficient**, taking values greater than or equal to 1 [@gomes_individual_2022;@neipel_powerlaw_2020] (see Appendix Equation \ref{power.trans}). The classic model with no heterogeneity and a linear transmission function corresponds to $\lambda=1$. Derivations of key mathematical results for SEIR models with heterogeneous mixing are given in the Appendix. For given values of the basic reproduction number ${\cal R}_0$ and $\lambda$, the final size of an unmitigated epidemic can be calculated from Appendix Equation \ref{eqn:fsizehet}. # Results ## 1. Predicted size of an unmitigated epidemic Figure \ref{fig:esize} shows that for any given value of ${\cal R}_0$, the predicted final size of an unmitigated epidemic is smaller with a model that allows for heterogeneity than with a model with no heterogeneity. Figure \ref{fig:tf} compares the predicted trajectory of the epidemic under a model with no heterogeneity ($\lambda=1$) with the trajectory under a model with heterogeneity quantified by $\lambda=4$, assuming ${\cal R}_0=2.5$. With $\lambda=4$, the herd immunity threshold is `r round(100 * (1 - 1 / 2.5^(1/4)))`%, and the epidemic ends when `r round(100 * finalsize(2.5, 4))`% have been infected. For comparison, the values for the size of an unmitigated epidemic predicted by the IC modellers (who assumed ${\cal R}_0 = 2.4$) and the LSHTM modellers (who assumed ${\cal R}_0=2.7$) were respectively 81% and 85%. These values are slightly smaller than the calculated value for the size of an epidemic under a model with no heterogeneity, because the IC and LSHTM models allowed for slight heterogeneity attributable to variation of contact rates between age bands, equivalent to using $\lambda=1.2$ in Appendix Equation \ref{eqn:fsizehet}. The underlying principle can be described without mathematics. Early in the epidemic, the most susceptible and the most highly connected individuals are selected for infection. Individuals who have persistently high contact rates with others -- the hubs of the social network -- are not only more likely to acquire infection but also more likely to transmit to others. As these highly susceptible or highly connected individuals are removed from the susceptible compartment, the infection will not spread so easily because the remaining individuals are more resistant to infection or have lower contact rates. The epidemic eventually kills itself by immunizing the most susceptible and the most highly connected individuals -- the hubs of the social network. Varying infectiousness (for instance where a few "superspreaders" infect many people) does not affect the final size of the epidemic unless infectiousness is correlated with susceptibility. ### Comparison of model-based predictions with the observed trajectory of the first wave In most countries the proportion infected in the first wave up to August 2020 was less than 10% [@vaselli_seroprevalence_2021; @rostami_sarscov2_2021]. Even in Iran, the first country outside China to be severely affected, seroprevalence was only 17% at the end of April 2020 [@poustchi_sarscov2_2021]. As contact rates fell in all these countries, there is no example of an unmitigated first wave. Seroprevalence rates of more than 50% in the first wave have been reported only from the Amazonian region: 77% in blood donors in Manaus [@buss_threequarters_2021] in October 2020 and 70% in Iquitos in July 2020 [@alvarez-antonio_seroprevalence_2021]. The Manaus blood donors result appears to be an outlier, far higher than the adjusted seroprevalence of 29% in a cohort of Manuaus residents studied in October 2020 [@buss_threequarters_2021]. ### Quantifying heterogeneity We can distinguish two sources of heterogeneity: biological susceptibility, and connectivity. There is indirect evidence for heterogeneity of biological susceptibility before the pandemic, based on studies showing that 20-50% of unexposed individuals had cross-reactive T cell responses to SARS-CoV-2 antigens [@sette_preexisting_2020]. Although later studies suggested that these cross-reactive T cell responses protect against infection [@swadling_preexisting_2022;@kundu_crossreactive_2022], no study has directly quantified the relation of infection risk to T cell reactivity in samples taken before the pandemic. In principle, the contribution of heterogeneity of connectivity to the immunity coefficient $\lambda$ could be calculated from population-based surveys of the variance of contact rates between individuals. However as surveys of contact rates have recorded only a single day's contacts in each individual [@mossong_social_2008;@gimma_changes_2022], no estimate of variance in persistent connectivity can be made since we can not separate between and within individual variability. The variance between age bands can be estimated however, and contributes an immunity coefficient of about 1.2 [@gomes_individual_2022]. If biological susceptibility contributes to heterogeneity, we expect new epidemic waves to occur when new variants arise that overcome pre-existing resistance to the dominant strain. If connectivity contributes to heterogeneity, we expect new waves to occur when social networks are rewired, for instance at the beginning of a new school year. We can distinguish six epidemic waves in Britain up to early 2022: wave 1 starting in March 2020, wave 2 in September 2020, wave 3 in December 2020, wave 4 in May 2021, wave 5 in September 2021, and wave 6 in December 2021. Waves 2 and 5 were not obviously related to new variants, and coincided with the start of the school year. Waves 3, 4, and 6 were attributable to new variants: respectively Alpha, Delta and Omicron. This suggests that both sources of heterogeneity contributed to the wave-like pattern of the epidemic. ### Learning the transmission function Under variability in susceptibility or mixing rate models there is a one-to-one correspondence between the distribution of susceptibility or mixing and the transmission function. Hence we do not need to estimate the distribution of susceptibility or mixing rates if we can learn the transmission function from the trajectory of the epidemic. The reproduction number $\cal R$ can be estimated from the rate of growth of the epidemic if the distribution of the serial interval -- the time between successive cases in a chain of transmission -- is known, without any other modelling assumptions [@anderson_reproduction_2020]. In the UK the reproduction number is estimated to have fallen from about 3 at the start of the epidemic to about 0.5 in mid-April 2020, staying at around 0.8 until late summer [@wood_inferring_2021;@wood_was_2021]. It is estimated that contact rates fell by 60% to 80% [@vaselli_seroprevalence_2021;@liu_rapid_2021;@gimma_changes_2022] after lockdown in comparison with the pre-pandemic period in the UK and other European countries. Two studies have attempted to estimate immunity coefficients from the trajectory of the first wave, with models that allowed for changing contact rates [@gomes_individual_2022; @tkachenko_timedependent_2021]. Using survey data to adjust for contact rates, the immunity coefficient $\lambda$ was estimated as 2.9 for Scotland and England [@gomes_individual_2022]. Using mobility data to adjust for contact rates, estimates of $\lambda$ ranged from 4.1 to 4.7 in cities and states of the USA [@tkachenko_timedependent_2021]. Although these results suggest that heterogeneity of connectivity is far more than is allowed for in the SAGE models, the results must be treated with considerable caution. It is statistically difficult to distinguish the effects of falling contact rates, reduced susceptibility among those remaining uninfected, seasonal effects on transmissibility, and the substantial modification of contact networks caused by lockdowns, because all these factors were changing at the same time. That said, @gomes_individual_2022 includes an extensive review of contact studies from which contact distribution parameters could be directly estimated: these suggest that the estimates were realistic. Although published reports from SAGE examined the sensitivity of their forecasts to various modelling assumptions, sensitivity to the assumption of a linear transmission function was not included in these analyses. The rationale for ignoring unmeasured heterogeneity is not clear: in comments on social media, modellers have expressed scepticism that heterogeneity can be modelled without measuring susceptibility directly or relying on strong assumptions about the distribution of susceptibility in the population. However where modelling results are highly sensitive to a process, the difficulty of measuring that process is a poor excuse for ignoring it. In any case, it is not necessary to make strong assumptions about the distribution of susceptibility if we can learn the transmission function from the observed trajectory of the epidemic, although it is clear that attempting to do so may well increase the level of uncertainty that has to be acknowledged in the modelling results. Conversely, if no reliable estimate of the transmission function can be made from real-world observations, that would seriously undermine the case for using models to set policy at all. ## 2. Requirement for critical care. The MRC Biostatistics Unit estimates that about 10% of the English population was infected by the end of the first wave [@birrell_realtime_2021]. Over the same period about 0.1% of the population died from COVID-19 as underlying cause, giving an infection fatality rate of about 1%. This is close to the value of 0.9% assumed in the IC model at least for the first wave, though the infection fatality rate was heavily weighted by the high fatality rate in care home residents, who accounted for half of all fatal cases in Scotland [@mckeigue_rapid_2020]. However, the SAGE modellers' assumptions about the ratio of cases requiring critical care to fatal cases were wide of the mark. The IC report assumed that the ratio of cases entering critical care (30% of 4.4% hospitalised) to fatal cases (0.9% of infections) would be 1.3 [@ferguson_report_2020]. From the LSHTM report [@davies_effects_2020] the value assumed for the ratio of cases entering critical care to fatal cases at the peak of the epidemic can be calculated (from the projections of peak requirement of 200,000 critical care beds with average length of stay 10 days, and peak weekly deaths 57000) as 2.5. In Scotland during the first wave the actual ratio of cases entering critical care to fatal cases was 0.2 [@mckeigue_rapid_2020]. As critical care capacity was not exceeded, this low ratio of cases entering critical care to death was presumably because most of those who were destined to die from COVID-19 were assessed as unlikely to benefit from critical care because they were very frail and near the end of life. Although the inability of health-care systems to cope with the first wave of the epidemic in Lombardy was interpreted as the result of failing to control population-wide transmission, a subsequent report noted that in Lombardy at this time "SARS-CoV-2 became largely a nosocomial infection" [@boccia_what_2020]. Nosocomial transmission accounts for a high proportion of severe cases because even if the clinically vulnerable can shield themselves from other sources of exposure, they cannot easily avoid exposure to hospital when they need medical care [@mckeigue_relation_2021a]. From 2 April 2020 onwards SAGE recorded increasing concern about nosocomial infection, but on 7 May noted that "Granular data are not yet currently available from PHE to fully understand transmission pathways in healthcare settings" [@scientific_advisory_group_for_emergencies_scientific_2020]. ## 3. Suppression or mitigation? The optimal mitigation scenario considered by the IC modellers -- 75% effective shielding of the 15% of the UK population who were aged over 70, combined with case isolation -- was predicted to reduce mortality by two-thirds. The modellers argued that this would not be enough to keep the requirement for critical care within the limits of capacity, and that this would require reducing the reproduction number below 1 to suppress the epidemic. Possibilities for a broader mitigation strategy based on risk-stratified protection of the vulnerable were not examined at this time. Based on the information available at the time about risk stratification, if all those past retirement age or with designated risk conditions had been advised to shield, at least 25% of the population would have been advised to shield initially, but this 25% at highest risk would have accounted for at least 85% of those at risk of fatal disease [@mckeigue_evaluation_2020]. Those who are advised that they are at high risk are likely to restrict their contact rates voluntarily: focused protection of this group can be supported without coercion. For those in the high-risk group who were economically inactive, not living with economically active adults and able to live without personal care, only minimal support for shielding would have been required. Additional support for those who needed it -- for instance because of exposure at work, carer responsibilities, or sharing a household with other adults who were likely to be exposed -- could have been provided at far lower cost to society than population-wide restrictions on social and economic activity. When the NHS Volunteer Responders Programme was announced on 24 March 2020 with the aim of supporting 2.5 million individuals identified as clinically vulnerable to COVID-19, 750,000 volunteers came forward within the first six days before recruitment was paused [@royal_voluntary_service_findings_2020]. By end of September 2020, though 385,000 volunteers had been approved and had used an app to register themselves as on duty, only 110,000 vulnerable individuals had been provided with support through this initiative. ## 4. Did suppression require lockdown? The modellers used survey data on contact rates to predict how non-pharmaceutical interventions would change the reproduction number. If transmissibility and susceptibility are constant, the reproduction number is proportional to the average contact rate, or to the largest eigenvalue of the matrix of age-stratified contact rates [@anderson_reproduction_2020]. The only survey of contact rates in the UK available in early 2020 was the POLYMOD survey in which 1012 individuals in the UK had recorded their contacts for a single day in 2005-6 [@mossong_social_2008]. As POLYMOD included only seven participants aged over 75, data on mixing in the most vulnerable age groups was wholly inadequate. The type of contact most relevant to aerosol transmission of respiratory viruses like SARS-CoV-2 -- sharing a poorly ventilated indoor space [@asadi_aerosol_2019] -- was not recorded: contacts were recorded only as physical (skin-to-skin contact) or non-physical (two-way conversation only). On 14 April 2020 SAGE noted that "Evidence suggests that transmission risk outdoors is significantly lower than indoors" but this did not lead to any reassessment of forecasts based on survey data that did not distinguish indoor contacts from outdoor contacts [@scientific_advisory_group_for_emergencies_scientific_2020]. As contacts were recorded only on a single day, the variance of persistent connectivity between individuals could not be estimated. It is doubtful whether any reliable predictions of the effect of interventions on the reproduction number could be made from such inadequate data. The attribution of the fall in reproduction number below 1 between March and April 2020 in the UK to lockdown [@flaxman_estimating_2020;@knock_key_2021] was based on very strong modelling assumptions: in a reanalysis with relaxation of these assumptions, the reproduction number was estimated to be falling rapidly before lockdown was imposed [@wood_inferring_2021;@wood_was_2021]. This was supported by reconstruction of infection dates from symptom onset dates reported in the REACT-2 survey. A similar analysis showed that the reproduction number had been falling before the second lockdown in November 2020, and this was supported by direct estimates of infection rates from the ONS COVID-19 infection survey [@wood_was_2021]. For Sweden, one of the few countries that did not impose a lockdown, the IC modellers predicted on 26 March that given ${\cal R}_0 = 2.7$, mitigation with "social distancing of whole population" without a lockdown would lead to 34895 deaths from COVID-19 [@walker_report_2020]. In the year up to the end of July 2020, the number of COVID-19 related deaths reported in Sweden was 5741 and the estimated number of excess deaths from all causes was 4329 [@rypdal_estimation_2021]. These results, together with studies from other countries casting doubt upon the necessity of lockdown were available well before the second lockdown in the UK (early May 2020 for @wood_inferring_2021) though final peer-reviewed publication was later [@kuchenhoff_analysis_2021;@chin_effect_2021;@bendavid_assessing_2021]. # Discussion The analyses above show that of the four propositions on which the recommendation for lockdown was based, one -- the assumption that 2% of those infected would require critical care -- was unequivocally wrong, and another -- that mitigation through focused protection would not be effective in limiting morbidity and mortality -- was not seriously questioned. The other two propositions -- that in an unmitigated epidemic 80% would be infected, and that only a lockdown could suppress the epidemic -- were reliant on strong but unrealistic modelling assumptions and weak data. Unlike weather models, epidemic models were not developed as forecasting models, but rather as theoretical models to aid the understanding of epidemic processes. While weather models are continuously validated and re-trained against measurements (and related climate models have been subject to extensive refinement of component processes via ground-truth measurement and back-cast validation), such a firm grounding in empirical reality is elusive in epidemic modelling and impractical for models of a newly emergent disease. At the same time, unlike weather and climate models, epidemic models are not based on well-understood physics, readily susceptible to accurate mathematical description, but rather on complex human behaviour, which at best can be captured only crudely by tractable mathematical models. In this respect epidemic models perhaps better resemble economic forecasting models, but without the substantial streams of relevant data available to tune the latter. It follows that epidemic models may well be suitable for answering such broad scale qualitative questions as 'if we coercively suppress all social contact to the maximum extent possible will we drive ${\cal R}<1$?', but are simply not designed or calibrated for more nuanced quantitative questions, such as 'how limited a set of restrictions would have a reasonable chance of avoiding severe health system overload?'. Given these limitations, although it may be appropriate to use epidemic models for exploration of what would happen if reality resembled the model, this exercise needs to be accompanied by proper communication of uncertainties and a greater willingness to rapidly re-assess the models when data suggest poor calibration. For example, given the strong sensitivity of epidemic wave size to variance in mixing rates and susceptibility, no prediction or statement of uncertainty should be considered reliable if it fails to take this variability into account. It also follows that studies aimed at *measuring* such variability should be a priority if model predictions are to be used in a policy setting. For example, a cross-sectional or cohort study with regular testing for infection, like the ONS COVID-19 Infection Survey [@pouwels_community_2021], could be used to quantify the heterogeneity of connectivity, if cohort members recorded diaries of contacts of different types. Similarly, if peripheral blood mononuclear cells were stored at baseline, such a cohort study could also be used to quantify the heterogeneity of susceptibility attributable to pre-existing T cell cross-reactivity. The modelling reports at the outset of the epidemic in March 2020 failed to communicate that by relying on the unrealistic assumption of no unmeasured heterogeneity they were likely to overestimate the size of the epidemic, and that no reliable prediction of the effects of non-pharmaceutical interventions could be made without more information about the mode of transmission. This suggests that policy advice in future epidemics should rely less on *models*, with a greater priority given to the rapid establishment of high quality direct *measurement* studies, such as the ONS survey, and on shoe-leather epidemiology [@koo_snow_2010] investigating the key target variable: transmission to vulnerable individuals. SAGE minutes from April-May 2020 show that although they were alerted to evidence that transmission risk was higher indoors than outdoors, and that nosocomial transmission was increasing, this did not lead to any serious reconsideration of policies. A report from NERVTAG and EMG to SAGE on 22 July 2020 concluded that "it is possible that transmission through aerosols could happen where a person who generates significant amounts of virus is in a poorly ventilated space with others for a significant amount of time." but gave only cursory attention to investigations of outbreaks that by this time had provided compelling evidence for aerosol transmission [@scientific_advisory_group_for_emergencies_nervtag_2020]. The choice between risk-stratified protection of the vulnerable and intervention to suppress population-wide transmission is not binary: clinically vulnerable individuals will need support to shield themselves whether or not attempts are made to suppress population-wide transmission. Both lockdowns and focused protection strategies divide the population into compartments with differing connectivity. With a lockdown, the high-connectivity compartment comprises essential workers and carers, and those who share a household with them. With focused protection of the vulnerable, the high-connectivity compartment would comprise those who are young, fit and not sharing a household with someone vulnerable. Sustainable measures to reduce transmission may also be preferable to "circuit-breaker" lockdowns [@keeling_precautionary_2021] in that they avoid repeated disruption of social networks that may lead to network rewiring with new hubs. Although herd immunity induced by natural infection of highly susceptible or highly connected individuals is likely to be only transient, even this would gain time for development of vaccines or other measures to protect the more vulnerable. Finally, the epidemic models treat all negative effects of interventions as externalities. This is not rational when the aim of epidemic management is to minimise all loss of life from the epidemic and associated disruption. It is also an example, along with ignoring heterogeneity, seasonality, nosocomial transmission and the aspects of transmissibility not captured by contact rates, of excluding effects that are difficult to quantify, irrespective of their importance for relevant outcomes. The consequent danger is of one-sided caution in the advice offered, where recommendations that are presented as precautionary with respect to the outcomes considered in the model, may be recklessly risky in terms of outcomes that are omitted. For example, although the Bank of England has estimated that the economic shock to the UK from the response to COVID is the largest in 300 years, and the effect of economic deprivation and economic shocks on mortality are well established [@marmot_health_2020], quantifying the effects of lockdown on years of life lost is not amenable to simple mathematical models of the kind used by the SAGE modellers. Given that the primary issue here is one of balance and fairness between controlling COVID and other health problems, the standard cost-effectiveness analyses usually employed to help ensure health equity offer an alternative approach. Although these still require some estimate of the life years saved by lockdowns (and other measures) simple bounding calculations are possible. These suggest that the cost per year of life gained by lockdowns was far higher than the threshold of \pounds30,000 used by NICE to decide whether to recommend use of a new pharmaceutical intervention in the health service [@miles_stay_2021]. # Declarations ## Availability of data and code This paper uses only publicly available data. R code to generate the figures and numbers in the manuscript is provided in the Markdown source document. ## Declaration of conflicting interests The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. ## Funding The authors received no financial support for this work. # References
\newpage ```{r episize, echo=FALSE, message=FALSE, fig.height=3, fig.width=5,fig.align='center', fig.cap="\\label{fig:esize}Relation of final size of unmitigated epidemic to immunity coefficient, for different values of the basic reproduction number $R$. See Appendix eqns. (2) and (3)."} esize <- function(lambda=1,R0=2.5) { if (lambda!=1) f <- function(x) (1-x)^(1-lambda) - (1 +x*(lambda-1)*R0) else f <- function(x) 1 - exp(-R0*x) - x #x0 <- optimize(f,c(0,1),maximum=TRUE)$maximum x0 <- .01 uniroot(f,c(x0,1))$root } es <- Vectorize(esize,"lambda") lambda <- seq(1,5,length=100) R0 <- c(2, 3, 4, 5) esize <- NULL for (i in 1:length(R0)) { esize.r0 <- data.table(lambda=lambda, fsize=es(lambda, R0[i])) esize.r0[, R0 := R0[i]] esize <- rbind(esize, esize.r0) } esize[, R0 := factor(R0)] setorder(esize, R0) ggplot(esize, aes(x=lambda, y=fsize, color=R0, linetype=R0)) + geom_line() + guides(color=guide_legend(reverse = TRUE)) + guides(linetype=guide_legend(reverse = TRUE)) + labs(color = expression(italic("R")[0]), linetype = expression(italic("R")[0])) + xlab("Immunity coefficient") + ylab("Final proportion infected") + scale_y_continuous(expand=c(0, 0), limits=c(0, 1)) ``` ```{r plotbyS, echo=FALSE, message=FALSE, fig.asp=1, fig.width=5, fig.cap="\\label{fig:tf}Relation of reproduction number and daily cases to the fraction of the population that is still susceptible, for immunity coefficients of one (no heterogeneity) and four. Vertical dotted lines are at one minus the herd immunity thresholds."} p.RS <- ggplot(data=epidem.dt[uncorr==TRUE], aes(x=St, y=R, color=lambda)) + geom_line() + theme(legend.position=c(0.8, 0.7)) + labs(color="Immunity coefficient") + scale_x_reverse(limits=c(1, 0), expand=c(0, 0)) + xlab("Susceptible fraction of population") + ylab(expression(paste("Reproduction number", italic("R")))) + scale_y_continuous(limits=c(0, 3.2), expand=c(0, 0)) + geom_hline(yintercept=1, linetype="dotted") + geom_segment(data=epidem.H, aes(x=H.S, xend=H.S, y=0, yend=1, color=lambda), linetype="dotted", size=1) p.dSS <- ggplot(epidem.dt, aes(x=St, y=cases, color=lambda)) + geom_line() + scale_x_reverse(limits=c(1, 0), expand=c(0, 0)) + theme(legend.position="none") + labs(color="Immunity coefficient") + xlab("Susceptible fraction of population") + ylab("Daily cases as fraction of population") + geom_vline(data=epidem.H, aes(xintercept=H.S, color=lambda), linetype="dotted", size=1) gA <- ggplotGrob(p.RS) gB <- ggplotGrob(p.dSS) grid::grid.newpage() grid::grid.draw(rbind(gA, gB)) #pdf("plotbyS.pdf") #grid::grid.draw(rbind(gA, gB)) #dev.off() ``` \newpage ```{r plotbydays, echo=FALSE, message=FALSE, fig.asp=1.5, fig.width=5, fig.cap="\\label{fig:plotbydays}Relation of daily cases to time since start of epidemic at different levels of heterogeneity, given basic reproduction number R_0=2.5. Immunity coefficient of 1 is equivalent to no heterogeneity. Dotted lines are herd immunity thresholds"} p.days <- ggplot(epidem.dt, aes(x=t, y=cases, color=lambda)) + geom_line() + theme(legend.position=c(0.8, 0.7)) + labs(color="Immunity coefficient") + scale_x_continuous(limits=c(0, 100)) + xlab("Days since prevalence reached 1/1000") + ylab("Cases as daily fraction of population") + geom_vline(data=epidem.H, aes(xintercept=H.t, color=lambda), linetype="dotted", size=1) p.R.days <- ggplot(epidem.dt, aes(x=t, y=R, color=lambda)) + geom_line() + theme(legend.position=c(0.8, 0.7)) + labs(color="Immunity coefficient") + scale_x_continuous(limits=c(0, 100)) + xlab("Days since prevalence reached 1/1000") + ylab("Reproduction number") + geom_vline(data=epidem.H, aes(xintercept=H.t, color=lambda), linetype="dotted", size=1) p.logcases <- ggplot(epidem.dt, aes(x=t, y=log10(cases), color=lambda)) + geom_line() + theme(legend.position=c(0.8, 0.8)) + labs(color="Immunity coefficient") + scale_x_continuous(limits=c(0, 100)) + scale_y_continuous(limits=c(-5, -1)) + xlab("Days since prevalence reached 1/1000") + ylab("Daily cases (log scale)") + geom_vline(data=epidem.H, aes(xintercept=H.t, color=lambda), linetype="dotted", size=1) gA <- ggplotGrob(p.days) gB <- ggplotGrob(p.R.days) gC <- ggplotGrob(p.logcases) #grid::grid.newpage() #grid::grid.draw(rbind(gA, gB, gC)) p.logcases.S <- ggplot(epidem.dt, aes(x=St, y=log10(cases), color=lambda)) + geom_line() + theme(legend.position="none") + labs(color="Immunity coefficient") + scale_x_reverse(limits=c(1, 0), expand=c(0, 0)) + scale_y_continuous(limits=c(-5, -1)) + xlab("Susceptible fraction of population") + ylab("Daily cases (log scale)") + geom_vline(data=epidem.H, aes(xintercept=H.S, color=lambda), linetype="dotted", size=1) #p.logcases.S ``` ```{r plotfsize, echo=FALSE, fig.width=5, fig.cap="Final size of epidemic at R0=3"} p.fsize <- ggplot(data=fsize.het, aes(x=lambda, y=finalsize, colour=factor(R0))) + geom_line() + scale_y_continuous(limits=c(0, 1), expand=c(0, 0)) + xlab("Immunity coefficient") + ylab("Final proportion infected") + theme(legend.position=c(0.6, 0.75)) #p.fsize ``` # Appendix: S(E)IR models with heterogeneity of susceptibility or mixing The key results are concisely derived here. For more see [@novozhilov_spread_2008;@gomes_individual_2022;@neipel_powerlaw_2020;@tkachenko_timedependent_2021]. ## Epidemic size under heterogeneity in susceptibility versus linear transmission (no heterogeneity) An S(E)IR model that allows susceptibility to vary between individuals can be written \begin{align*} \frac{d ~}{d t} s(\alpha,t) &= - \alpha s(\alpha,t)I(t)\\ &\cdots \\ \frac{d I}{d t} & = \cdots -\gamma I(t) \end{align*} where $\cdots$ indicates details of model structure that do not change the results here. The positive transmission rate parameter $\alpha$ varies over the population (with susceptibility). Note that other authors write the transmission rate as the product of a population-level parameter $\beta$ representing transmissibility, and an individual-level susceptibility parameter $\omega$ defined to have initial mean of 1 [@novozhilov_spread_2008]. $s(\alpha, t) d\alpha$ is the fraction of the population who are in the susceptible compartment at time $t$ and have parameter value betwen $\alpha$ and $\alpha + d\alpha$. So if $S_t$ is the fraction of the population who are susceptible at time $t$ then $s(\alpha, t)/S_t$ is the probability density function of $\alpha$ within the susceptible compartment. Without loss of generality we define $S_0$ to be 1. The first equation can immediately be integrated $$ \int_{s(\alpha,0)}^{s(\alpha,t)} \frac{1}{s} ds = - \alpha \int_0^t I(t^\prime) dt^\prime \Rightarrow \log \left \{\frac{s(\alpha,t)}{s(\alpha,0)} \right \} = -\alpha q_t \text{ where } q_t = \int_0^t I(t^\prime) dt^\prime. $$ Hence $s(\alpha,t) = s(\alpha,0)\exp(-\alpha q_t)$. $q_t$ has value 0 at $t=0$, and increases with $t$. Now suppose that the initial probability distribution for $\alpha$ is $\text{gamma}(k,\nu)$ with p.d.f. $s(\alpha,0) = \nu^{k}\alpha^{k-1}e^{-\alpha \nu}/\Gamma(k)$. It follows directly that $s(\alpha,t) = \nu^{k}\alpha^{k-1}e^{-\alpha (\nu + q_t)}/\Gamma(k)$. To obtain $S_t$ we need to integrate $s(\alpha, t)$ over $\alpha$ (from $0$ to $\infty$ where unspecified), $$ S_t = \int s(\alpha,t)d\alpha= \frac{\nu^k}{\Gamma(k)} \int \alpha^{k-1}e^{-\alpha (\nu + q_t)} d\alpha = \frac{\nu^k }{(\nu + q_t)^{k}}. $$ The second integral uses the identity \[ \int \frac{(\nu + q_t)^{k} \alpha^{k - 1} e^{-\alpha (\nu + q_t)}}{\Gamma(k)} d\alpha \equiv 1, \] since the p.d.f. of a $\text{gamma}(k,\nu + q_t)$ distribution integrates to one (like any p.d.f.). Similarly to find the time derivative of $S_t$, we must integrate $\alpha$ out of the expression for the time derivative of $s(\alpha,t)$. $$ \frac{dS}{dt} = -I(t) \int \alpha s(\alpha,t) d \alpha = -I(t) \frac{ \nu^k}{\Gamma(k)} \int \alpha^k e^{-\alpha (\nu + q_t)} d\alpha = - I(t) \frac{ \nu^k }{\Gamma(k)} \frac{\Gamma(k + 1)}{(\nu + q_t)^{k+1}}. $$ Again the second integral is obtained from the fact that a $\text{gamma}(k+1,\nu + q_t)$ p.d.f. must integrate to one. Using the identity $k\Gamma(k) \equiv \Gamma(k+1)$ and defining $\lambda = 1 + 1/k$, routine manipulation gives the transmission function as: \begin{equation} \frac{dS}{dt} = - \frac{k}{\nu }S_t^{1+1/k} I(t) = - \frac{k}{\nu} S_t^\lambda I(t). \label{power.trans} \end{equation} ${\cal R}_0$, the expected number of new infections produced by an infectious individual at the start of the epidemic, is the mean $\alpha$ at time zero multiplied by the expected duration of infectiousness. The mean of a $\text{gamma}(k,\nu)$ random variable is $k/\nu$ so ${\cal R}_0 = k / (\nu \gamma)$. The final epidemic size is found by dividing (\ref{power.trans}) by the rate of removal $dR/dt = \gamma I$, to get $dS/dR = - {\cal R}_0 S_t^\lambda$. Rearrangement and integration gives $$ \int_1^{S^{\infty}} S^{-\lambda} dS = - \int_0^{R_\infty} \mathcal{R}_0 dR \Rightarrow \frac{S_\infty^{1-\lambda} - 1}{1 - \lambda} = -R_\infty \mathcal{R}_0 $$ The final proportion $x$ infected is $R_\infty = 1 - S_\infty$, so routine re-arrangement then shows that $x$ must satisfy \begin{equation} x = 1 - [ 1 + (\lambda - 1) \mathcal{R}_0 x ]^{-1 / ( \lambda - 1 ) } \label{eqn:fsizehet} \end{equation} For $\lambda=1$ the same approach shows that $x$ must satisfy \begin{equation} x = 1 - \exp{\left( -\mathcal{R}_0 x \right)} \label{eqn:fsizehom} \end{equation} This is the limiting case in which $k \to \infty$ and $\nu \to \infty$, while $k/\nu \to$ a constant (for any meaningful model). For this to happen the variance of $\alpha$ must tend to zero - corresponding to no variability in susceptibility. ## Individual variation in mixing rates The preceding model assumes that infectiousness is independent of susceptibility. This is a poor model for the case in which variability in transmissibility is related to variability in social contact rates. Therefore consider a model in which each individual has their own value of parameter $\alpha$, and the transmission probability between individuals with parameter values $\alpha$ and $\alpha^\prime$ is proportional to $\alpha\alpha^\prime$. Defining $I(\alpha,t)$ as the equivalent to $s(\alpha, t)$ for the infectious class, the susceptibility model is then modified to \begin{align*} \frac{d ~}{d t} s(\alpha,t) &= - \int \alpha s(\alpha,t) \alpha^\prime I(\alpha^\prime,t) d\alpha^\prime\\ & = -\alpha s(\alpha,t) \bar \alpha_t^\prime I(t) \end{align*} where $\bar \alpha_t^\prime = \int \alpha^\prime I(\alpha^\prime,t)/I(t) d \alpha^\prime$. To proceed analytically, assume that the infectious state is short enough that, to a good approximation, the distribution of $\alpha$ in the infectious state at time $t$ is given by the distributions of $\alpha$ in those first becoming infectious at time $t$. That is $I(\alpha,t) \propto \alpha s(\alpha,t)$. Hence $\bar \alpha_t^\prime = \int \alpha^2 s(\alpha,t) d\alpha/\int \alpha s(\alpha,t) d \alpha$ (the dividing integral is the normalizing constant for the p.d.f. of $\alpha$ in the $I$ state). Now redefining $q_t = \int_0^t \bar \alpha^\prime_{t^\prime} I(t^\prime) dt^\prime$, the algebra follows through exactly as in the variable susceptibility case so that $s(\alpha,t)=s(\alpha,0)e^{-\alpha q_t}$, and again assuming a $\text{gamma}(k,\nu)$ initial $\alpha$ distribution we get $$ \frac{dS}{dt} = - \frac{k}{\nu} S^{1+1/k} \bar \alpha_t^\prime I(t). $$ The need to evaluate $\bar \alpha_t^\prime$ is new, but straightforward given that we again have an explicit gamma-like form for $s(\alpha,t)$. Using the same gamma integral tricks as above, $$ \bar \alpha_t^\prime = \frac{\Gamma(k+2)(\nu + q_t)^{k + 1}}{\Gamma(k+1)(\nu + q_t)^{k+2}} = S_t^{1/k} \frac{k+1}{\nu}. $$ Clearly $\bar \alpha_0^\prime = (k+1)/\nu$, and the product of this, the initial expected value of $\alpha$, and the expected residence time in the $I$ compartment gives the basic reproductive number ${\cal R}_0 = k(k+1)/(\nu^2\gamma)$. Hence \begin{align*} \frac{dS}{dt} &= - \frac{k(k+1)}{\nu^2} S_t^{1+2/k} I(t) \\ &= - {\cal R}_0 \gamma S_t^\lambda I(t) \end{align*} where now $\lambda = 1 + 2/k$. Clearly the equation to solve for the final epidemic size is again (\ref{eqn:fsizehet}), all that has changed is the relationship between $\lambda$ and $k$. ## Other distributions for susceptibility and contact rates It is easy to work with any initial $\alpha$ distribution for which the moment generating function $\mathcal{M}(s) = \mathbb{E}(e^{s \alpha})$ exists (a very mild restriction). For a wide range of distributions the m.g.f., its inverse and derivatives are known and readily computed. Under either model we have \begin{equation} S_t = \int s(\alpha,0) \exp(-\alpha q_t) d\alpha = \mathcal{M}(-q_t), \label{mgf1} \end{equation} as the integral is simply the expectation defining the m.g.f. Differentiating the integral w.r.t. $q_t$, and writing $\mathcal{M}^\prime(\cdot)$ for the first derivative of $\mathcal{M}$, we also have, $$ \mathcal{M}^\prime(-q_t) = \int \alpha s(\alpha, 0) \exp(-\alpha q_t) d\alpha $$ Hence, for the variable susceptibility model the transmission function is $$ \frac{dS}{dt} = - \int \alpha s(\alpha, t) d\alpha I(t) = -\mathcal{M}^\prime(-q_t) I(t) = -\mathcal{M}^\prime\{\mathcal{M}^{-1}(S_t)\}I(t) $$ where $\mathcal{M}^{-1}$ is the inverse function of $\mathcal{M}$ and $-q_t = \mathcal{M}^{-1}(S_t)$ follows directly from (\ref{mgf1}). So again we have a one-dimensional ODE. $h(S) = \mathcal{M}^\prime\{ \mathcal{M}^{-1}(S)\}$ is known as the *effective susceptible fraction*. The variable mixing model is equally straightforward in the general distribution case. All that changes is that the right hand side of the ODE for $S_t$ is again multiplied by the term $\bar \alpha^\prime_t$. But differentiating the defining equation for $\mathcal{M}(-q_t)$ once more and substituting into the definition of $\bar \alpha^\prime_t$ we obtain $\bar \alpha^\prime = \mathcal{M}^{\prime\prime}(-q_t)/\mathcal{M}^\prime(-q_t)$, and hence the transmission function is $$ \frac{dS}{dt} = - \mathcal{M}^{\prime\prime}\{\mathcal{M}^{-1}(S_t)\} I(t), $$ another one dimensional ODE with $h(S) = \mathcal{M}^{\prime\prime}\{\mathcal{M}^{-1}(S_t)\}$. @novozhilov_spread_2008 gives more detail, including examples for several distributions, and argues that the power transmission function obtained by assuming an initial gamma model is a good approximation more widely. However, once the original models are reduced from infinite dimensional to systems of a few readily computed ODEs in this way, they can readily be solved numerically to any desired accuracy, so approximation is somewhat superfluous.