What the reproductive number R0 can and cannot tell us about COVID-19 dynamics.

The reproductive number R (or R0, the initial reproductive number in an immune-naïve population) has long been successfully used to predict the likelihood of pathogen invasion, to gauge the potential severity of an epidemic, and to set policy around interventions. However, often ignored complexities have generated confusion around use of the metric. This is particularly apparent with the emergent pandemic virus SARS-CoV-2, the causative agent of COVID-19. We address some misconceptions about the predictive ability of the reproductive number, focusing on how it changes over time, varies over space, and relates to epidemic size by referencing the mathematical definition of R and examples from the current pandemic. We hope that a better appreciation of the uses, nuances, and limitations of R and R0 facilitates a better understanding of epidemic spread, epidemic severity, and the effects of interventions in the context of SARS-CoV-2.

Introduction

With the emergence of SARS-CoV-2, the novel coronavirus responsible for COVID-19, much attention has been given to the reproductive number, , and its initial state, (Viceconte and Petrosillo, 2020). is the expected number of infections generated by an infected individual in an otherwise fully susceptible population and in the absence of interventions (Anderson and May, 1991, Diekmann et al., 1990). Since each infection produces an average of new infections, describes the exponential growth of infections during the early phase of an epidemic. Under relatively general assumptions, can be used to determine the probability an emerging disease will cause an epidemic, the final size of an epidemic, and what level of vaccination would be required to achieve herd immunity (Anderson and May, 1991, Delamater et al., 2019, Heffernan et al., 2005, Roberts, 2007). In many cases, has been invaluable such as for predicting the risk of measles resurgence (Béraud et al., 2018, Hens et al., 2015) and for managing emerging infectious diseases such as SARS-CoV-1 (Lipsitch et al., 2003) and foot-and-mouth disease (Ferguson et al., 2001). When interpreted correctly, and in conjunction with additional relevant information, it can yield valuable insight. However, misinterpretation may lead to faulty conclusions regarding disease dynamics.

The virus SARS-CoV-2 emerged in Wuhan, China in late 2019 and has since become pandemic causing over one million deaths worldwide by September 30th, 2020 (Johns Hopkins University & Medicine Coronavirus Resource Center, 2020) in addition to severe economic distress. Policy makers have relied on estimates of to tailor control measures (e.g. Ferguson et al., 2020), but these estimates vary tremendously within and between populations around the globe (Fig. 1). It is important to understand why these estimates vary. It is also important to understand how the utility of the reproductive number is limited. Here, we derive and explain some of the key nuances of and , paying particular attention to insights and limitations with respect to the emerging pathogen SARS-CoV-2.

Fig. 1

Estimates of the of SARS-CoV-2 vary substantially between locations (Binny et al., 2020, Chen et al., 2020, Choi and Ki, 2020, Deb and Majumdar, 2020, Giordano et al., 2020, Johndrow et al., 2020, Ke et al., 2020, Korolev, 2021, Lewnard et al., 2020, Li et al., 2020b, Majumder and Mandl, 2020, Mizumoto et al., 2020, Peirlinck et al., 2020, Pitzer et al., 2020, Ranjan, 2020, Read et al., 2020, Riou and Althaus, 2020, Sanche et al., 2020, Senapati et al., 2020, Shim et al., 2020, Singh and Adhikari, 2020, Tang et al., 2020, Wu et al., 2020, Xiao et al., 2020, Yuan et al., 2020, Zhao et al., 2020). Each point represents a literature-compiled average estimate for a different geographic area (sample size noted alongside means, error bars show plus or minus 1 standard error). For individual studies that provided multiple estimates for a single geographic area, the median estimate was used to avoid pseudo-replication. An analysis (not shown) confirmed that estimation method (transmission model, exponential growth model, or stochastic simulation method) did not drive the pattern of variation in by location. Recent meta-analyses of values for SARS-CoV-2 consider the effects of estimation methods in more detail (Alimohamadi et al., 2020, Barber et al., 2020).

Deriving and

A general definition

How many new infections will be caused by a single infected individual? For a directly transmitted pathogen, the answer to this question can be written as: Above, is the reproductive number, is the rate of contacts that an infected individual has with susceptible individuals at time post infection, is the probability that a contact at time results in a new infection, and is the probability of still being infected at time . Notably, we could have combined and into a single parameter since the probability of infection given contact falls to zero after an individual recovers, but we prefer this more explicit formulation. Eq. (1) yields , the total number of infections one infected individual would generate over the course of their infection. When the population is fully susceptible and when no preventative interventions have been imposed, as would be expected at the beginning of a novel outbreak, Eq. (1) yields . Note that we have neglected to explicitly incorporate individual variation and temporal variation in contact rates, the probability that a contact results in a new infection, and the time to recovery. However, and , are intended to be averages, and so this variation is inherently a part of the reproductive number calculation. This variation could be explicitly included with additional subscripts to denote all possible infected hosts, noninfected hosts, times since the epidemic began, and ages of infections. For simplicity, we will not use any subscripts for , , and when we discuss how they relate to the reproductive number of SARS-CoV-2 in the remaining sections of this commentary.

Calculating the reproductive number using an epidemiological model

Estimating the individual parameters in Eq. (1) requires extensive data collection for a specific pathogen, host population, and time. Alternative approaches to estimating and therefore frequently rely on epidemiological models and epidemic data.

Following the lead of Kermack and McKendrick (1927), epidemics are often modeled as a set of ordinary differential equations such as for the classical Susceptible–Infected–Recovered (SIR) model. To illustrate the derivation of the reproductive number, we present a model of SARS-CoV-2 transmission that captures some of the key characteristics of the system. Unlike in the classic SIR model, in our model, infected individuals can be asymptomatic or not, and all non-asymptomatic infections progress through two stages, a pre-symptomatic and a symptomatic stage. Transmission is possible from each of these infected classes, but at different rates. Individuals showing symptoms are detected through testing (e.g. fever checks or diagnostic testing), and these individuals are quarantined. We assume no transmission from the quarantined class, and we assume that only symptomatic or quarantined individuals die from infection. Otherwise, all infected individuals eventually recover. In our model, the state variables , , , , , and are respectively the densities of susceptible, asymptomatically infected, pre-symptomatic, symptomatic (subscript “” used to denote COVID-19), quarantined, and recovered individuals in a population. Note that “” here is distinct from the reproductive number “”, but we use both for historical reasons.

Here, the subscripts , , , and denote parameters for asymptomatic, pre-symptomatic, COVID-19 symptomatic, and quarantined individuals respectively. The parameters are transmission coefficients, the parameters are virus-induced death rates, and the parameters are recovery rates. is the fraction of new infections that are asymptomatic, is the infection progression rate from pre-symptomatic to symptomatic, and is the detection rate of symptomatic infections. Parameter values can, in principle, be estimated from data (e.g. Ferretti et al., 2020). Note that further complexity could be added to this model, for example, to incorporate spatial structure, hospitalization, age structure, host-specific exposure risk, superspreading, or time-varying mortality. Simulations of our model are shown in Fig. 2.

Fig. 2

Simulation of transmission model of SARS-CoV-2 presented in Eqs. (2.1)–(2.6) for (left) and (right) . We achieved this change in by altering the parameter , which is the per day detection rate of symptomatic cases, thus illustrating the potential value of rapid case detection. Note that a larger leads to a faster epidemic and a greater fraction of the population becoming infected. The horizontal, dotted gray line in each panel indicates the herd immunity threshold (i.e. the value of that corresponds to ). Note that infection declines, but does not immediately disappear, after crossing the herd immunity threshold. Parameter values: ,(left) , (right) . Initial conditions: .

Numerous methods can be used to derive from a model (Heffernan et al., 2005). Perhaps most famously, is the dominant eigenvalue of the next generation matrix (Diekmann et al., 1990). In the following text, we use the survival function to calculate because of its intuitive connection to Eq. (1), but if we had used the next generation matrix, we would have derived the exact same . The survival function uses three components to determine the number of new infections caused by an initial case: (1) the rate at which an individual in a particular class causes new infections (from Eq. (1), ), (2) the probability that an individual is still in the class at time (from Eq. (1), ), and (3) the probability that an initial infected individual will enter that class. Note that this third term is implicitly assumed to be one in Eq. (1). The integral of the product of the first two terms multiplied by the third term yields the contribution to that comes from a focal class. Deriving the overall requires summing over all possible classes. To illustrate, consider the asymptomatic contribution for an average initial infection. Here, (1) is equal to , (2) is equal to , and (3) the probability that the initial case enters the asymptomatic class is equal to . Using these values, and extending the same methods to all infectious classes we derive: where is the time since entering class . Since changes slowly at the beginning of an epidemic when few individuals are infected, we can treat as a constant with respect to time. This assumption allows us to analytically solve Eq. (3.1), which yields Evaluating Eq. (3.2) at time rather than time 0 yields . Notice that the reproductive number has contributions from the asymptomatic, pre-symptomatic, and symptomatic classes, but not from the other classes because the other classes cannot cause new infections. If transmission were possible from these classes, additional terms would need to be included in Eq. (3.2). Alternatively, if we simplified our model by assuming that all infections were symptomatic () and pre-symptomatic individuals were not infectious (), then Eq. (3.2) would simplify to .

Calculating the reproductive number without an epidemiological model

can also be estimated without an epidemiological model, which can be especially useful if parameter estimates or even an appropriate model structure is not yet known. In principle, one could calculate by simply counting the cases attributed to infected individuals at or near the beginning of an outbreak. In practice, this method is rarely employed since contact tracing networks are rarely established during the earliest phase of an emerging disease outbreak (but see Pung et al., 2020) and estimates could be inaccurate due to bias towards observing large chains of transmission.

can also be inferred from the growth rate of cases early in an outbreak. Since the number of susceptible individuals changes slowly during the initial stages of an outbreak, early case growth rates can be approximated by exponential growth: the number of cases , where is the epidemic growth rate. If the number of cases is known for at least two time points, one could calculate the epidemic growth rate . The relationship between and depends on the distribution of the generation interval , which is defined as the amount of time between infection of two individuals where the second infection is caused by the first. can be approximated by direct observation or specified by an epidemiological model (Wallinga and Lipsitch, 2007). For the classical SIR model, the generation interval is exponentially distributed, leading to (Wallinga and Lipsitch, 2007), where is the mean of . Other relationships can be calculated for other distributions of (Wallinga and Lipsitch, 2007). See Zhao et al. (2020) for an example of this method used to calculate the of SARS-CoV-2. It is important to recognize however that emerging epidemics often grow more slowly than exponential due to stochastic effects, small population sizes, network effects, or preventative measures (Chowell et al., 2016). In these cases, a generalized growth model may more accurately reflect epidemic growth (Viboud et al., 2016).

After epidemics have started (i.e. when populations are no longer fully susceptible and when interventions may have been imposed), it is still useful to calculate to understand if cases will continue to grow or decline. A recent review explored various methods for this calculation (Gostic et al., 2020). Gostic et al. (2020) recommended two methods. The first, developed by Cori et al. (2013) uses time series incidence data and the serial interval to calculate in real time, which is particularly useful for assessing the impact of interventions as they are employed. The second, developed by Wallinga and Teunis (2004) uses similar data and is useful for the retrospective calculation of .

No matter the method used to calculate it, limited data or unreliable data early in an epidemic can make it difficult to constrain . The World Health Organization originally estimated the of SARS-CoV-2 to be between 1.4 and 2.5 (WHO, 2020). More recent estimates of have varied from 2.2 to 6.47 for the beginning of the Wuhan outbreak (Fig. 1). This represents tremendous uncertainty when attempting to use for public health planning. For example, if we were using these estimates to design a vaccine campaign capable of achieving herd immunity for a vaccine with perfect efficacy, our vaccination target (calculated as under assumptions of Eq. (2.1)–(2.6)) would be 29% of the population at or 85% at . Even with better estimates of , however, misconceptions around this metric lessen its practical utility.

Misconception 1: The reproductive number alone explains future dynamics

As we have explained, is calculated during the early stages of an epidemic because of its value in determining future infection dynamics in the absence of intervention (Anderson and May, 1991, Ma and Earn, 2006). is similarly useful for assessing changes in transmission over time and the potential impacts of interventions. However, and cannot fully explain future dynamics under particular circumstances. We discuss two important situations in which their value is limited: first, when hosts become aware of infection and alter their behavior, and second, when individual infection risk is heterogeneous.

Shifts in behavior that influence contact rates or the probability of infection given contact can alter over extremely short timescales (see Eq. (1)). For example, as awareness of the SARS-CoV-2 epidemic grew in the United States in March 2020, human mobility ground to a near halt (Gao et al., 2020, Warren and Skillman, 2020), presumably reducing contact rates . Other individual behavioral changes such as increased handwashing and mask wearing (Belot et al., 2020, Goldberg et al., 2020) have reduced the probability of transmission given contact (Eikenberry et al., 2020, Liang et al., 2020). Such bottom-up forces combined with top-down government-imposed interventions (e.g. school closures, banned gatherings) reduced to below 1 (the threshold for epidemic persistence) by late April 2020 in some states (Johndrow et al., 2020, Miller et al., 2020). Similar reductions to were documented in China (Li et al., 2020a, Tian et al., 2020) and other countries (Ensser et al., 2020, Giordano et al., 2020, Kupferschmidt, 2020, Yuan et al., 2020). Indeed, in models of the 1918 influenza pandemic, incorporating a behavioral response to death rates improved model fits (Bootsma and Ferguson, 2007, He et al., 2013). For SARS-CoV-2, behavioral changes have caused to fluctuate above and below 1 at different times based on the perceived threat of COVID-19 (Santamaría and Hortal, 2021).

While behavioral changes can temporarily reduce as described above, more sustainable reductions in are typically achieved when susceptible individuals are removed from populations either through naturally acquired immunity or vaccination. When hosts are heterogeneous such that some individuals are more likely to contract infection than others, declines faster than predicted by Eq. (2.1)–(2.6) or other SIR models lacking host heterogeneity (May and Anderson, 1987). This is because those most susceptible (for example, due to high exposure or low inherent immunity) will become infected earlier in an epidemic, leaving individuals that are on average more resistant (Langwig et al., 2017, May and Anderson, 1987). Thus, the classical formulation of the “herd immunity” threshold, , does not apply to populations with heterogeneity in risk of infection. To illustrate this point, Britton et al. (2020) developed a COVID-19 model with heterogeneous mixing between age classes and additional heterogeneity in contact rates. Under their assumptions, herd immunity was reached after 43% of the population acquired immunity through natural exposure, rather than the 60% that would be required under their assumptions in the absence of heterogeneity. Other studies have predicted that heterogeneity could have an even larger effect if susceptibility was correlated with infectivity (Gomes et al., 2020, Tkachenko et al., 2020). In some of the hardest hit areas, an appreciable fraction of people have been infected with SARS-CoV-2 (e.g. 22.7% in New York City, USA (Rosenberg et al., 2020) and over 44% in Manaus, Brazil (Buss et al., 2020)). However, the impact of these levels of infection on future disease dynamics is unknown, since for SARS-CoV-2, heterogeneity in infection risk is still highly uncertain (Randolph and Barreiro, 2020). Moreover, heterogeneity is likely to change through time due to changes in individual or government-mandated responses (Dolbeault and Turinici, 2020). Thus, even after herd immunity has been reached, additional waves of infection could occur if heterogeneity changes (Tkachenko et al., 2020).

Misconception 2: The reproductive number is constant over space

The of many pathogens are often referred to as known values. For example, the of measles is 12–14, polio is 5–7, and pertussis is 12–17 (Doherty et al., 2016). For SARS-CoV-2, estimates typically range from 2–3 (Liu et al., 2020a). However, the parameters (, , and in Eq. (1)) that make up can differ substantially from place to place (Figure 1, Delamater et al., 2019). It follows that interventions to reduce to less than 1 may need to vary in aggressiveness across locations (Stier et al., 2020).

Since differs between groups of people, combining multiple groups together to estimate a population-wide can produce misleading notions of disease spread. For example, a high average measles vaccination rate in the United States keeps below 1 nation-wide, but localized communities with high rates of vaccine refusal still experience serious outbreaks (Leslie et al., 2018). In the current COVID-19 pandemic, disease transmission has so far been much higher in refugee and low income populations compared to non-refugee and high income populations (Chopra and Sobel, 2020, Lau et al., 2020, Ruiz-Euler et al., 2020) and has been documented to vary between US counties (Sly et al., 2020). Since is an average, combining communities with high and low transmission may yield an estimate of , yet disease may still readily spread (Li et al., 2011). On the other hand, splitting populations may mean missing transmission events that occur between populations, thus underestimating (Smith et al., 2009).

Awareness of the consequences of how people are grouped can help us interpret values of the reproductive number. Within groups, behavior associated with work, home, and recreation affects contact rates . As we have discussed above, SIR models often assume that contact rates, and thus and , depend on host density. Although, evidence is mixed as to whether larger cities have higher values of and for SARS-CoV-2 (Heroy, 2020, Stier et al., 2020), built environments (e.g. hospitals, airport terminals, factories) do often have high values of (Dietz et al., 2020). This may partially explain patterns of explosive SARS-CoV-2 transmission in venues such as cruise ships and meat packing facilities (Althouse et al., 2020, Dyal et al., 2020, Mizumoto and Chowell, 2020). Household contacts have also been particularly important to the transmission dynamics of the virus (Bi et al., 2020). Therefore, differences in between populations can be partially explained by differences in household sizes between countries and cultures (Singh and Adhikari, 2020). For example, large, multigenerational households in Italy may have contributed to the large and deadly outbreak there relative to other European countries where multigenerational households are less common (Dowd et al., 2020). Similarly, the probability a new infection results from contact and the probability of remaining infected over time may vary by population. For SARS-CoV-2, individuals with severe symptoms have 60 times more viral RNA in nasal swabs, which likely increases both their ability to transmit the virus and the amount of time they remain infected (Liu et al., 2020b). Since older individuals are more likely to develop severe infection (Yang et al., 2020), is likely to be greater in populations with older individuals, such as in nursing homes (McMichael et al., 2020) or in developed countries (Dowd et al., 2020).

Misconception 3: The reproductive number is enough to tell us how large an epidemic will be

It is tantalizing to imagine that can be used to predict the extent of an outbreak, since it can be calculated during the early stages of an epidemic. Indeed, is related to final epidemic size (Kermack and McKendrick, 1927), but this relationship can be substantially affected by the fraction of the population infected initially and by heterogeneity in transmission.

If we rescale population sizes such that the initial susceptible population size , and we assume that population sizes are sufficiently large to neglect demographic stochasticity (Hartfield and Alizon, 2013, Tildesley and Keeling, 2009), then the fraction of the population infected during an epidemic can be determined from the final epidemic size equation. This equation, , is not sensitive to the epidemic model used to calculate so long as the model assumes that the population is well-mixed with homogeneous susceptibility. Note that this equation prominently features , the fraction of the population infected at the beginning of the outbreak (or at the beginning of an intervention). Efforts to reduce the reproductive number below 1 are understandably a high priority, but when is close to or less than 1, the final outbreak size is more sensitive to changes in the fraction of infected individuals than it is to changes in (Fig. 3). This is because more infected individuals will fuel the outbreak for longer, infecting a greater proportion of the susceptible population even when each case produces on average less than one new infection. Now that SARS-CoV-2 infection rates are already substantial in populations around the globe, must be considered in addition to . For example, Pei et al. (2020) estimated that implementing social distancing policies one week earlier could have reduced the cases in the United States by early May, 2020 by 55% (over 700,000 cases) by keeping the number of infected individuals low at the time such policies were implemented.

Fig. 3

Epidemic size contours and shading show that when is close to 1, the epidemic is more strongly influenced by a reduction of than by a reduction of . For instance, if (dashed line), the epidemic could infect from less than 0.1% to greater than 16% of the population as ranges from just above 0% to 3% of the population. Epidemic size was calculated using the final size equation, , where . Shading indicates the cube root of epidemic size with lighter colors corresponding to smaller outbreaks.

While a final epidemic size can be calculated using and , the final size equation above does not apply to populations with heterogeneous infection risk (Andreasen, 2011, Ball, 1985, Hébert-Dufresne et al., 2020, Ma and Earn, 2006). Given the inherent heterogeneity in human social networks, it is surprising that the final size equation is so accurate for many diseases such as childhood diseases (Caudron et al., 2015). This may be because transmission networks among children are close to random (Bjørnstad et al., 2002) or have other features that allow the final size equation to be accurate (Bansal et al., 2007). For diseases where this is not the case (e.g. many sexually transmitted diseases), heterogeneity could in principle be incorporated into the final epidemic size equation (Dwyer et al., 2000), but estimating heterogeneity early in an epidemic can be challenging. Moreover, as we describe in misconception 1, heterogeneity in infection risk can change over time as a result of human behavior or interventions, such as the shutdowns in response to the COVID-19 epidemic (Dolbeault and Turinici, 2020, Ruiz-Euler et al., 2020). Estimates of how future government restrictions and behavioral changes will alter heterogeneity in infection risk are thus critical for assessing the likely impact of the outbreak (Gomes et al., 2020). Such estimates are also key in determining thresholds for herd immunity (Randolph and Barreiro, 2020) and in prioritizing the distribution of interventions such as vaccines when they first become available (Atkinson and Cheyne, 1994, Giambi et al., 2019).

Conclusions

As we have discussed, the reproductive number and its initial value can be used to assess the potential for disease invasion and persistence, to predict the extent of an epidemic, and to infer the impact of interventions and of relaxing control measures. Determining in real time is especially helpful for the latter goals. For example, the of SARS-CoV-2 was used to justify implementations of lockdowns in the United Kingdom (Ferguson et al., 2020), and is one factor that is considered in lockdown relaxation policies (Thompson et al., 2020).

However, the utility of and can easily be overstated. Though we have focused on three misconceptions that we felt were particularly important for the COVID-19 pandemic, there are additional considerations to appreciate for accurate interpretation of and (Heffernan et al., 2005, Li et al., 2011, Roberts, 2007). These include complications associated with stochasticity (Keeling and Grenfell, 2000), superspreading (Lloyd-Smith et al., 2005), metapopulation dynamics (Cross et al., 2007), seasonality (Bjørnstad et al., 2002), multiple hosts (Roberts, 2007), and pathogen evolution (Hartfield and Alizon, 2014). Some of these nuances may also be important for SARS-CoV-2 dynamics. For example, superspreading events appear to be common with this virus (Althouse et al., 2020). Although superspreading affects epidemic dynamics, masks this feature of disease spread because it only represents average transmission.

The reproductive number is a valuable tool for understanding disease dynamics. describes early disease growth, and it can be estimated when little else is known about a pathogen. Later, can be used to measure the effectiveness of interventions. However, a singular focus on and can lead to inaccurate conclusions. Considering the nuances of disease dynamics when interpreting the reproductive number allows for a stronger understanding of the patterns with which the SARS-CoV-2 virus has traversed the globe, why it has impacted some populations more than others, and how best to limit future transmission.