Acceleration of plague outbreaks in the second pandemic.

Historical records reveal the temporal patterns of a sequence of plague epidemics in London, United Kingdom, from the 14th to 17th centuries. Analysis of these records shows that later epidemics spread significantly faster ("accelerated"). Between the Black Death of 1348 and the later epidemics that culminated with the Great Plague of 1665, we estimate that the epidemic growth rate increased fourfold. Currently available data do not provide enough information to infer the mode of plague transmission in any given epidemic; nevertheless, order-of-magnitude estimates of epidemic parameters suggest that the observed slow growth rates in the 14th century are inconsistent with direct (pneumonic) transmission. We discuss the potential roles of demographic and ecological factors, such as climate change or human or rat population density, in driving the observed acceleration.

Plague epidemics have afflicted human populations since at least the sixth century (1, 2). These events have had dramatic and long-lasting effects on human demography and behavior, especially those outbreaks associated with the second pandemic (14th to 19th centuries) in Europe and Asia (1, 3–5), and have inspired many theoretical studies of the ecology and evolution of infectious disease (6–13). We are now in the third pandemic (Modern Plague), with outbreaks continuing to occur in some parts of the world (14–18). Plague also remains a source of concern due to the bioterror potential of the causative agent, Yersinia pestis (19, 20).

Recent advances in paleogenomics have definitively established that historical plague pandemics were caused by Y. pestis (21, 22), as proposed in the 19th century after Yersin discovered the bacterium’s link to bubonic plague (23). Researchers have reconstructed the evolutionary history of plague and other pathogens by sequencing and reconstructing nearly complete pathogen genomes from persistent DNA fragments (21, 24, 25). The strain isolated from victims of the Black Death (London 1348) is remarkably similar to extant human strains (Modern Plague): the core genomes*of these strains are similar (21), which makes it challenging to identify important evolutionary or ecological patterns from genomic investigations alone. Here we complement genetic studies by exploring more traditional (historical, demographic, and epidemiological) sources of information from a 300-y span of plague outbreaks in the same location (London), revealing a striking change in plague transmission dynamics over the course of the Renaissance period, namely, a fourfold increase in the initial growth rate of outbreaks.

We quantify this change without making any assumptions about the underlying transmission processes, exploiting methodology that we have developed previously for this purpose (26). We then consider how this inference can contribute to the debate concerning whether plague transmission was primarily indirect (via rat fleas) or direct (pneumonic human-to- human). We argue that strictly pneumonic transmission in the 14th century is implausible but that beyond this the best that can be done at present is to highlight the biological complexities and uncertainties that limit the potential for further inferences.

Data

The city of London, United Kingdom, is unusual in the extent to which patterns of death and disease can be reconstructed from extant documents. We have analyzed three sources of data (Fig. 1).

Fig. 1.

(Top Left) Part of a will proved in the PCC, dated 18 December 1644 (34). (Top Center) A parish register page from St Giles without Cripplegate, August 1665 (38). Image credit: Wellcome Collection, licensed under CC BY 4.0. (Top Right) One of the LBoM, for the week beginning 26 September 1665 (photo by Claire Lees, taken in the Guildhall Library, City of London). (Bottom) Mortality in London, United Kingdom, 1340 to 1380 and 1540 to 1680, aggregated 4-weekly, plotted on a log scale. The three distinct sources of data () are last wills and testaments of Londoners whose wills were probated in the Court of Husting (14th century) or the PCC (16th to 17th centuries), weekly aggregations of burials listed in extant parish registers (29), and weekly plague deaths listed in the LBoM (27). Major plague years are highlighted in yellow. Aside from these and various minor plague years, there are notably unusual patterns during an epidemic of sweating sickness in 1551 (ref. 29, p. 70), the influenza epidemic of 1557 to 1559 (ref. 29, p. 70) (which coincided with the end of the reign of Bloody Mary I and the ascension of Elizabeth I in 1558 [indicated by a crown icon]), and the absence of a monarch during the Interregnum from 30 January 1649 to 29 May 1660.

London Bills of Mortality.

The London Bills of Mortality (LBoM) (27, 28) provide weekly counts of deaths by cause during each of the known plague epidemics in the 16th and 17th centuries. The LBoM thus give us information specifically about mortality attributable to plague for each epidemic, which we will use as a proxy for plague incidence (6, 26). When available, the LBoM are our most reliable records (except for the epidemic of 1593; ).

Parish Registers.

After 1538 (ref. 29, p. 54, and ref. 30), London parish registers provide information on mortality. We use weekly counts of deaths published online by Cummins, which are estimated to cover 80% of total deaths (31). During plague epidemics, most deaths were likely due to plague; we also adjust for nonplague mortality by estimating a baseline death rate from the years before and after an epidemic ().

Wills and Testaments.

Before 1538, no direct tabulations of mortality for London are available. However, we do have detailed information on last wills and testaments, which can be used as proxies for mortality (29). In particular, we have digitized and tabulated daily counts of wills recorded in the Court of Husting (32) for the 14th century; these give us data for the plague epidemics of 1348, 1361, 1368, and 1375 [monthly counts for three of these epidemics were previously published by Cohn (33)]. We use online data from the Prerogative Court of Canterbury (PCC) (34) for the plague epidemics of the 16th and 17th century. Plague epidemics in London are mentioned in historical documents from the 15th century (ref. 27, chap. IV), but during this period, wills recorded in the Courts of Husting and Canterbury were too sparse to enable identification of epidemics.

Two dates are associated with each will: the date on which it was written and the date on which it was probated (i.e., accepted by the court after the death of the testator). Traditionally, demographers have used annual counts of wills by probate date as testators are known to have died by these dates. While this approach works well for inferring long-term trends in annual mortality, probate dates are inappropriate for reconstructing detailed time series of mortality over short time scales as there is a variable and irregular lag between the (unknown) date of the death and the probate date; most of the wills associated with 14th century epidemics were probated several months after they were written (35) ().

Instead, we use the counts of wills written in a given time interval as a proxy for plague incidence. The date on which a will was written may precede the testator’s death by a long period, but during severe epidemics, will dates were likely to have been correlated with the fear of infection (and hence with disease incidence), just as internet searches for influenza symptoms can predict 21st century epidemic patterns (36).

When using wills as proxies for mortality or incidence, we must keep in mind that 1) a much smaller proportion of the population is sampled (the will counts in Fig. 1 are about 10 times smaller than the counts based on parish registers) and 2) the subpopulation of individuals who wrote wills is strongly biased toward merchants and others who owned property (primarily males).

For the period since 1540, the existence of multiple sources allows us to cross-check results—in particular, to test our assertion that aggregated counts of wills by date written are an adequate proxy for mortality rates—and avoid relying on any questionable or poorly sampled segments of data (37) (). In 1563 and 1593, the number of PCC wills is too small to extract a signal from the noise, but the data for the later three outbreaks verify that the numbers of wills written follow the plague epidemic patterns evident in the mortality records (Fig. 2). Based on cross-correlation analysis and differences between epidemic peak times, will counts appear to have lagged mortality recorded in parish registers and LBoM by approximately 3 to 5 wk; our inferences about epidemic growth rates appear to be robust to this lag ().

Fig. 2.

Observed time series (points) and phenomenological model fits (lines; ) that yield the estimated growth rates (listed in Table 1 and plotted in Fig. 3). The data sources () are described in the main text and the legend of Fig. 1. For visual comparison, wills were aggregated weekly to match the frequency of mortality observations (fits to wills were based on the original daily counts). Vertical dashed lines show 1 April, 1 July, and 1 October of each year. Weeks with zero counts are shown along the bottom edge of the graph (present in all 14th century epidemics, 1625, and 1665). , display the data during the major plague epidemics on linear rather than logarithmic scales.

Table 1.

Maximum likelihood estimates (MLEs) of the initial exponential growth rate () and doubling time () with their 95% CIs, obtained from the time series shown in Fig. 2 (see )

						Assuming pneumonic plague
Source	Epidemic	Growth rate(1/y)	Doubling time (days)	R_0,P	Attack rate Z_P	R^2
Husting wills	1348	4.5 (3.4, 6.2)	55.6 (41.1, 74.1)	1.06 (1.04, 1.09)	0.11 (0.08, 0.16)	0.439
Husting wills	1361	9.2 (4.9, 15.4)	27.5 (16.4, 51.6)	1.14 (1.07, 1.25)	0.23 (0.12, 0.37)	0.693
Husting wills	1368	8.4 (3.3, 21.7)	30.2 (11.7, 77.2)	1.12 (1.04, 1.38)	0.21 (0.08, 0.49)	0.333
Husting wills	1375	20.3 (6.0, 60.3)	12.5 (4.2, 42.2)	1.35 (1.08, 2.41)	0.47 (0.15, 0.88)	0.669
London bills	1563	21.8 (20.1, 23.7)	11.6 (10.7, 12.6)	1.38 (1.34, 1.42)	0.49 (0.46, 0.52)	0.961
London parish	1563	20.6 (18.9, 22.5)	12.3 (11.2, 13.4)	1.35 (1.32, 1.39)	0.47 (0.44, 0.51)	0.978
London parish	1593	15.9 (14.2, 17.3)	15.9 (14.6, 17.8)	1.26 (1.23, 1.29)	0.38 (0.35, 0.41)	0.986
Canterbury wills	1603	14.0 (8.3, 23.8)	18.0 (10.6, 30.4)	1.23 (1.12, 1.42)	0.34 (0.21, 0.53)	0.853
London bills	1603	29.1 (26.5, 31.8)	8.7 (8.0, 9.5)	1.54 (1.48, 1.60)	0.61 (0.57, 0.64)	0.927
London parish	1603	19.9 (18.2, 21.8)	12.7 (11.6, 13.9)	1.34 (1.30, 1.38)	0.46 (0.43, 0.49)	0.988
Canterbury wills	1625	27.1 (11.4, 53.5)	9.3 (4.7, 22.3)	1.49 (1.18, 2.19)	0.58 (0.28, 0.84)	0.900
London bills	1625	27.3 (25.9, 28.4)	9.3 (8.9, 9.8)	1.50 (1.46, 1.52)	0.58 (0.56, 0.60)	0.996
London parish	1625	23.0 (22.3, 23.7)	11.0 (10.7, 11.3)	1.40 (1.39, 1.42)	0.51 (0.50, 0.53)	0.999
Canterbury wills	1665	35.8 (17.6, 47.5)	7.1 (5.3, 14.4)	1.69 (1.29, 2.01)	0.69 (0.42, 0.80)	0.842
London bills	1665	22.2 (21.0, 23.4)	11.4 (10.8, 12.0)	1.39 (1.36, 1.41)	0.50 (0.48, 0.52)	0.984
London parish	1665	23.2 (20.6, 26.0)	10.9 (9.7, 12.3)	1.41 (1.35, 1.47)	0.52 (0.47, 0.56)	0.987
Husting wills	Early	5.9 (3.5, 9.8)	43.2 (25.7, 72.5)	1.08 (1.04, 1.15)	0.15 (0.08, 0.25)
Canterbury wills	Late	23.0 (11.7, 45.0)	11.0 (5.6, 21.5)	1.40 (1.18, 1.94)	0.51 (0.29, 0.78)
London bills	Late	23.9 (21.3, 26.7)	10.6 (9.5, 11.9)	1.42 (1.37, 1.48)	0.53 (0.49, 0.57)
London parish	Late	20.4 (18.3, 22.8)	12.4 (11.1, 13.8)	1.35 (1.31, 1.40)	0.47 (0.43, 0.51)

The goodness of fit measure () is the proportional reduction in the mean squared error, i.e., where is the model mean squared error and is the (population) variance of the data; predictions and observations are aggregated to weekly before computing the for wills, for consistency among data sources. The implied basic reproduction numbers () and attack rates (), assuming the generation interval distribution for Modern pneumonic plague (44), are also shown. All MLEs and associated CIs are shown in Fig. 3.

Fig. 3.

Estimated initial growth rates () and 95% profile CIs on a logarithmic scale, for each of the epidemics shown in Fig. 2. Second and third panels show mixed-model estimates () of overall average growth rate for each epoch (early, 14th century vs. late, 16th to 17th centuries) and data source. To aid interpretation, the additional vertical scales show the implied intrinsic reproductive number () and epidemic final size () under the assumption that the generation interval distribution (40) was the same as that estimated from 20th century pneumonic plague epidemics (44). All estimates and CIs are listed in Table 1 (which also lists the associated doubling times).

The LBoM indicate that there were 18 other plague epidemics in London between 1563 and 1666; these were all of much smaller magnitude and are discussed only in . No plague epidemics have occurred in London since the end of the Great Plague in 1666: a total of 77 deaths from plague are recorded in the LBoM after 1666, never more than 5 in a single week. The last plague deaths recorded in the LBoM occurred during the week of 6 May 1679.

Approach

Estimation of Epidemic Growth Rates.

While our collection of high-resolution mortality time series during plague epidemics allows us to estimate epidemic growth rates, we are constrained by the limitations of our data—the time series for each epidemic are short and noisy, especially during the 14th century (e.g., the epidemics of 1368 and 1375 peaked at a total of five wills written per week). Thus, we chose to use simple phenomenological models with a small number of parameters that can be fitted to short time series (starting size, growth rate, and total size; ) and do not attempt to separate different sources of variability or estimate the influence of various mechanistic processes [King et al.’s exhortation to use epidemic models that separate process from observation noise (39) applies primarily to forecasting, which we are not attempting]. Furthermore, we primarily consider the initial growth rate , which—unlike the basic reproduction number —can be estimated reliably from mortality data alone without any knowledge of the disease characteristics or natural history (26).

We estimate the epidemic growth rate for each epidemic as the maximum likelihood estimate of the initial growth rate of a logistic function fitted to the cumulative death curve (using first differences to avoid overconfidence; see for more detail), along with likelihood profile CIs. To summarize the overall difference between epochs, we then use the estimates for each epidemic as observations in a linear mixed model including fixed effects of epoch (early/14th century vs. late/16th to 17th centuries) and data source, and a random effect of outbreak year (see below), where each estimate was weighted by its certainty.

Finally, to robustly quantify the statistical significance of the difference between epochs, we ran a permutation test. We computed the null distribution of the between-epoch difference in growth rate by scrambling the estimates and associated CIs randomly across epochs. For each of the 126 possible permutations, we refitted the mixed model and extracted the estimated among-epoch difference in growth rate.

Estimating and Attack Rate.

Our analysis focuses on the epidemic growth rate , rather than the more commonly investigated basic reproduction number [the expected number of secondary cases caused by a primary case in a completely susceptible population (7)], because estimates of depend on information about the life history of the pathogen as well as on the epidemic curve itself (40). However, if we can estimate , then we can predict epidemic outcomes that we are unable to predict from alone; then, by comparing observed with theoretical outcomes under different scenarios, we can make inferences about modes of transmission.

In particular, if the host population is homogeneously mixed, we can use to estimate the expected final size (the proportion of the population infected by the end of an epidemic, also called the attack rate) (6, 41). The observed attack rate based on mortality records depends in turn on the fraction of infected individuals who die from disease. We can thus compare theoretical expectations of observed attack rates under different transmission modes with historical information about observed attack rates (see ).

We can use the distribution of the generation interval (the times elapsed between the moment a host is infected and the times at which they infect others) to compute from (40). If we know only the mean generation interval , we can still approximate , by assuming that the generation intervals are exponentially distributed; this approximation underestimates in the typical case where the generation interval distribution is less variable than an exponential distribution with the same mean.

Beyond the information on generation interval that we use to derive from , the observed attack rate depends on additional characteristics of an epidemic. If the population does not mix homogeneously, the final size will typically decrease (e.g., figure 6 in ref. 42). Since the case fatality proportion (CFP; the proportion of plague cases who die) is not 100%, the observed attack rate is less than the final size. Observed attack rates also decrease if some individuals are infected without showing symptoms (asymptomatic infections) (43) or without symptoms identified as plague (subclinical infections). Finally, the observed attack rate depends on the initial proportion of the population susceptible, which will be lower if some individuals have previously acquired immunity after prior nonfatal plague infections.

We can use independent estimates of the generation interval to estimate (and hence the theoretical attack rate) for London epidemics under different transmission modes. Given the limited information that we have, we consider the additional (extreme) assumptions that 1) the population was homogeneously mixed, 2) the population was initially 100% susceptible, 3) there were no asymptomatic infections, and 4) all infected individuals died (i.e., the CFP was 100%). All of these assumptions increase the observed attack rate (as detected from wills or deaths), so we can calculate an upper bound on the total mortality we would expect to observe in an epidemic with a particular transmission mode.

Results

The fitted models do as good a job as one could hope for in capturing the initial patterns of epidemic growth evident in the noisy will counts and agree very well with the patterns of growth in the mortality time series that are available in the later epoch (Fig. 2).

The growth rates for the 14th century epidemics are lower than those for the 16th to 17th centuries (Fig. 3); they may also be more variable than in the later period, although some of the appearance of variability is due to larger uncertainties (presumably due to the relatively poor sampling in the early epoch). The mixed model quantifies the dramatic increase in growth rate (“acceleration”) of epidemics between epochs. The growth rates increased fourfold in the late vs. early epochs: increased 3.9 [95% CI ] (). The permutation test results support the statistical significance of this difference; the observed order of outbreaks yields the most extreme difference in growth rates out of all 126 possible orderings, giving a two-tailed P value of .

Discussion

Plague epidemics in London appear to have been faster in the 16th to 17th centuries than in the 14th. Our analysis shows that this difference is very likely real, rather than driven by changes in the types of documentary evidence that are available; in particular, we have shown that inferences from will counts are consistent with those from death counts when both are available, in the 17th century.

Why the later plague epidemics were faster than the earlier ones is not clear. Nevertheless, to provide some context, we consider below a number of mechanisms that could have contributed to increases in epidemic growth rates. We then consider what our estimated epidemic growth rates suggest about the mode of transmission during the various London plague epidemics.

Possible Causes of Acceleration.

What factors could have caused plague epidemics in London in the 16th and 17th centuries to accelerate relative to the earlier plagues in the 14th century?

Pathogen or host evolution.

Faster growth could be a consequence of the pathogen evolving to infect hosts more effectively or to remain infectious for longer. However, such evolution may be challenging to demonstrate, particularly in light of the strong genetic similarity between the Second Pandemic and Modern Plague strains (21). Evolution of host resistance has also been suggested as a cause of changes in plague dynamics (11), although we have no particular mechanisms in mind for how the evolution of resistance over the course of centuries could accelerate epidemics.

Shifts in transmission mode.

Modes of transmission may have differed between epochs, which could account for differences in epidemic growth rates. The reservoir host of Y. pestis in London was rat populations where it was primarily spread by flea vectors such as Xenopsylla cheopis (ref. 14, p. 54). This transmission mode gives rise to bubonic plague in humans, as a consequence of zoonotic spillover; epidemic growth rate is driven by spread among rats and fleas, with human infection and mortality as a secondary consequence. In contrast, when Y. pestis infection in humans spreads to the lungs, it gives rise to pneumonic plague, which can then spread directly from human to human like other severe respiratory infections. It has also been more recently hypothesized that human ectoparasites (including lice and human fleas) could have supported indirect plague transmission without involving rats (13, 45).

Pneumonic plague probably has shorter generation intervals than bubonic plague (44, 46) and is often thought to spread much faster at the population level (44); one potential explanation for our observations is that 14th century epidemics were primarily bubonically transmitted, while 16th to 17th century epidemics were primarily pneumonic. With currently available data, we are unable to confirm or reject this hypothesis, but we consider below whether any inferences about transmission mode can be made at present.

Ecological and demographic changes.

The complete multispecies ecology of Y. pestis is extremely complicated (47). We focus on humans and not on the rest of the ecological community supporting Y. pestis because the only data available for the epochs studied here concern humans. Human population size and population density in London increased enormously from the 14th to 17th centuries (48–50) (), and density increases were exacerbated in the 17th century by the isolation of the sick and their contacts in pest houses (51). Since bubonic plague epidemics are driven by spillover from rat–flea epidemics, human population density would only have a direct effect on the speed of pneumonic epidemics. However, changes in human population (and living conditions) almost certainly affect rat, and flea, density (52). In addition, higher rat densities make it more likely that fleas departing dying rats end up on susceptible rat hosts. The magnitudes of these effects are extremely difficult to estimate, but in we explore them in a (crude) quantitative manner, without attempting to model the truly complex population dynamics that occur even during a single epidemic (46, 53).

Environmental change.

The external environment represents the third side of the host–parasite–environment disease triangle (54). Northern European climates changed significantly between the 14th and 16th to 17th centuries (the coldest period of the Little Ice Age occurred in the 17th century; ref. 55 and figure 1 of ref. 56). Furthermore, late-epoch epidemics occurred later in the calendar year and hence in different seasons (Fig. 2). A variety of studies link climate (e.g., overall wetness/dryness) or weather (e.g., annual precipitation or temperature) to plague incidence or transmission (57–63). Climatic changes may have been responsible for the changes in plague epidemics in London between 1348 and 1666, but it is challenging to make reliable inferences due to the lack of consensus about climate and weather variations in Medieval/Renaissance Europe (64).

Implications of Growth Rate Estimates for Transmission Mode.

One might hope to exploit ancient DNA (2, 21) to identify the transmission mode during London plague epidemics. Unfortunately, there are currently no identified genetic differences between strains that have caused indirect or direct transmission, so it will be difficult to determine from genetic analyses whether one mode or the other was dominant in any given Medieval/Renaissance epidemic.

Our growth rate estimates provide a starting point from which we can begin to explore the extents to which direct or indirect transmission played a role in the various London plague epidemics. Current data are too sparse and limited to make precise predictions (26, 65), but we consider the extremes of estimated parameter ranges that would favor each mode of transmission and ask whether the resulting predictions of attack rates are consistent with the rough knowledge of epidemic sizes available from the literature.

In the early epoch, primarily pneumonic transmission is implausible.

If we assume pneumonic transmission and use generation interval estimates from 20th century outbreaks (44) (i.e., assuming that the natural history of infection in London was similar to that of Modern pneumonic plague), in the early epoch we find [CI (1.04, 1.16)], which implies an attack rate of [CI (7.5%, 27%)] (left axes in Fig. 3). This estimate of total mortality is the largest possible for the early-epoch epidemics (violations of any of the four extreme assumptions discussed in would lower the total mortality). Total mortality of 15% is implausibly low in comparison to typical estimates for 14th century plague epidemics [e.g., “[b]etween the years 1346 and 1352, [plague] caused the death of … one-third of the world’s population at that time” (ref. 19, p. 163)] and certainly inconsistent with the careful and comprehensive analysis of Creighton for the 1348 to 1349 epidemic in London specifically [“the mortality was about one-half the population” (ref. 27, p. 128)]. It is therefore unlikely that 14th century plague epidemics in London were primarily pneumonic.

In the early epoch, growth rates appear to be consistent with bubonic plague.

What if we assume bubonic transmission? Unfortunately, the expected total human mortality from a bubonic plague epidemic depends on several poorly known aspects of rat ecology. For the generation interval for bubonic plague we use a rough estimate of 18 d ( y) based on the flea incubation period of 9 to 26 d (46), which dominates the time scales of other processes involved in rat–flea plague dynamics (see for further detail). In combination with a growth rate estimate of y, we obtain and . However, this estimate of the final size gives the fraction of rats, not humans, infected over the course of the epidemic. Translating this rat final size to the final size of the human epidemic requires us to quantify the amount of spillover from rats to humans. In particular, we need to know the rat-to-human ratio and the number of infected humans per infected rat, which can both be crudely estimated as 1:1 (53) (although these ratios will vary during the course of plague epidemics). Given an estimated CFP of 70 to 80% for early-epoch plagues (66), our (very rough) mortality predictions are reasonably consistent with the observed human mortality rates.

In the late epoch, primarily pneumonic transmission is unlikely but cannot be ruled out.

Now suppose that the late-epoch plagues were pneumonic with a generation interval similar to that of Modern pneumonic plague. The implied reproduction number (Fig. 3) would be [CI (1.2, 2.1)], which yields a maximum attack rate of 51% [CI (27%, 81%)], compared to the estimated of the population who died in each of the late-epoch epidemics (ref. 31, p. 4). Given that the CFP for (Modern) pneumonic plague is nearly 100% (20, 67–69), reconciling the predicted attack rates with observed mortality would require strong effects of heterogeneity, a large fraction of individuals resistant to plague due to prior nonfatal plague infections, a large proportion of asymptomatic or subclinical infections, or a decrease in transmission rate as the plague spread (perhaps due to behavioral avoidance of transmission). The effects of heterogeneity and behavioral avoidance may indeed have been stronger in the late epoch, as a result of concentration of the wealthy in central London and a tendency for them to flee London during plague epidemics (ref. 31, p. 4). While we still consider direct (pneumonic) transmission in late-epoch London plague epidemics unlikely because of the strong effects that would be needed to reduce the expected mortality rates to the observed levels, we cannot rule it out.

In the late epoch (as in the early epoch), growth rates appear to be consistent with bubonic plague.

Similar calculations with the estimated for the late epoch and the approximate generation interval for bubonic plague give estimated [CI ] and predicted rat attack rates of [CI ]. The human would be identical, if we were to make the same assumptions as above. This attack rate is even higher than the 51% predicted for pneumonic plague, but the CFP for bubonic plague was probably only 30 to 60% in the 16th to 17th centuries (e.g., refs. 19 and 69), which could bring the total mortality roughly in line with the observed . The lower CFP of bubonic plague would also lead to a build-up in the population of individuals who had survived previous plague epidemics and retained immunity (reducing the susceptible population available to be attacked), which in combination with heterogeneity and asymptomatic infections makes bubonic transmission plausible for the late epoch.

Other possibilities.

As mentioned above, it is possible that human ectoparasites could have supported indirect plague transmission (13, 45). While the biological details of this transmission mode are still uncertain, we can guess that the generation interval of this transmission mode would be intermediate between direct and indirect rat–flea transmission, giving rise to estimates of and intermediate between those discussed above. Investigations of plague biology in fleas have also raised the possibility that rat transmission may have shorter generation intervals than previously thought, which would lower estimated values for bubonic plague (70, 71) (but cf. ref. 72).

COVID-19 Context.

The world is currently in the midst of a pandemic of COVID-19 disease, caused by the SARS-CoV-2 virus (73). The infection fatality ratio (IFR) for COVID-19 is uncertain due to case detection limitations, but the evidence accumulated to date about this IFR (74, 75) indicates that it is substantially lower than for the 1918 influenza (76, 77) and much lower than for the historical plague epidemics studied here. Nevertheless, COVID-19 has had a much greater impact than seasonal influenza epidemics, which cause of order 500,000 deaths worldwide annually (76, 78).

Epidemic growth rates (and doubling times) have been important elements of the COVID-19 public health discourse (79). The quality and quantity of COVID-19 data far exceed that of the historical plague data we have analyzed in this paper. In particular, many countries report daily counts of newly reported cases, hospitalizations, and deaths. The techniques we have used to estimate initial growth of plague outbreaks are consequently easier to apply to modern data streams and should yield more accurate estimates. In addition, improved information about generation intervals will make it easier to convert growth rate estimates into estimates of the basic reproduction number . Of course, modern data are by no means perfect, and there are still many issues—including estimating generation intervals while accounting for disease dynamics, and assessing the importance of asymptomatic and presymptomatic transmission—that can lead to overconfident or incorrect estimates of disease parameters for COVID-19 (79–81).

Conclusion.

We have estimated epidemic growth rates for all of the recorded plague epidemics in London, United Kingdom, during the second pandemic (1348 to 1666). The major plagues of the 16th and 17th centuries grew much faster than the early plagues of the 14th century when the Black Death first invaded human populations. This conclusion is based solely on estimated growth rates—in particular, our estimates do not rely on assumptions about the mode of transmission or natural history of infection. The cause of this substantial acceleration is currently unclear.

We considered the implications of our estimated growth rates with respect to the mode of transmission and concluded that direct transmission (pneumonic plague) was unlikely to have been the primary mode in the early epoch. Definite conclusions about indirect transmission (bubonic plague) are difficult because of the substantial uncertainties about rat and flea ecology during this period.

We have created and assembled—and made available with this paper—machine-readable databases composed of thousands of handwritten records spanning more than 300 y. People living in London in 1348 or 1665 could not have imagined how these records might be used hundreds of years later. The time series that have emerged have allowed us to reveal and characterize large changes in patterns of infectious disease transmission. In concert with other forms of information (31, 82), these records should continue to contribute to a broader and more insightful view of population-level processes in early eras of history.

Materials and Methods

Growth Rate Estimates.

We have previously developed methods to use mortality data to estimate initial epidemic growth rates (26). The naïve approach of fitting a straight line to the relationship between logarithmic death counts and time at the beginning of the epidemic depends sensitively on the fitting window and produces overly narrow confidence intervals. Instead, we use maximum likelihood estimation to fit a curve that resembles the shape of a deterministic epidemic model (see ) and compute likelihood profile confidence intervals (83); this approach yields robust estimates of initial growth rates and accurate assessments of uncertainty. All variation around the fitted line is assumed to arise from observation errors. See Ma et al. (26) for methodological details.

In contrast to Ma et al. (26), who assumed a Poisson error model to fit simulations that were themselves generated with that assumption, the results shown here applied a negative binomial error model, which allows for the broader range of variation found in real data. We initially tried both logistic and Richards (84) models for the epidemic curve (i.e., the expected mean of the data) and both Poisson and negative binomial models for the distribution around the mean. These variants gave qualitatively similar conclusions about the pattern of growth rates, but we found that the Richards model was too unstable to fit reliably across all epidemics, so we reverted to fitting the logistic model for all epidemics. Similarly, the negative binomial model gave unstable results when the distributions were too close to Poisson; when the negative binomial dispersion parameter was estimated as , we reverted to a Poisson fit (see for further details of the optimization).

In order to summarize the differences between epochs, we fitted a linear mixed model [using the glmmTMB package in R (85)] to the estimates of log growth rates. To account for differing precision of the estimates across epidemics, we approximate the SE of the estimates as the width of the 95% confidence interval divided by 3.92 (the width in standard errors of a 95% normal CI), square it to find the variance, and scale the residual variance for each epidemic curve by this observed variance. This procedure is equivalent to weighting each epidemic curve by the inverse of its variance. We quantified the difference in epidemic growth rate between epochs based on the predicted difference in between wills in each epoch.

The permutation test P value (0.016) reinforces our confidence that late-epoch epidemics were faster than early-epoch epidemics; however, since it is larger than the parametric P values derived from the model fit (which are ), it suggests that the confidence intervals shown in Fig. 3 may be too narrow.

In light of recent work investigating the possibility of initial epidemic growth that is subexponential, we also considered fits based on the more general growth model of Chowell and coworkers (86–88). This approach leads to slightly slower initial growth estimates but does not change our conclusion that the later plagues were much faster than the earlier plagues.

Phenomenological Model.

Following Ma et al. (26), we model the cumulative mortality curve (or will-accumulation curve; henceforth, we refer to wills or death records as “mortality”) using a logistic model with an added baseline linear growth rate,Here is the initial cumulative mortality, is the initial exponential growth rate, and is the final size of the epidemic (). In the fitting procedure we use the unitless parameter ; all parameters are fitted on a transformed, unconstrained scale [i.e., , , and ]. The derivative is the instantaneous mortality rate. To avoid the inherent correlations between observations in the cumulative curve, we fit the observed data to the differences from Eq. 1: , where is the observation interval (1 d for wills and 1 wk for mortality data). The parameter represents baseline (nonplague) mortality. For the LBoM fits we set (since the LBoM data include only plague mortality); for the wills and parish register fits we estimated by computing the average observed mortality rate over the 2 y preceding and following the outbreak window (see ).

As discussed above, we initially fitted a Richards model, which extends the logistic equation with a shape parameter, but found that we could not get stable fits to the smallest/noisiest epidemics; we never found qualitative differences between growth rates estimated from logistic and Richards fits. Ma et al. (26) found that a delayed logistic curve, which adds an exponentially distributed delay between infection and mortality, estimated growth rates from mortality curves with less bias. However, the delay between infection and mortality in our data is extremely uncertain—it would represent another parameter that would have to be estimated from noisy data. For wills data, the delay could even be negative; i.e., wills could have been written out of fear of future infection. Furthermore, because we found little difference between the Richards and logistic fits—which Ma et al. (26) found to have opposite biases—growth rates estimated by these two methods are likely to bracket the growth rate estimated by the delayed logistic.

Fitting Windows.

How well a given model can fit a time series may depend on the time window used for fitting (88). In order to avoid choosing a window separately for each individual epidemic, we considered window selection criteria that could be applied uniformly to all epidemics without first examining the data. In all cases we defined the end of the fitting window to be one data point past the observed peak date (the date of maximum observed deaths). For wills and parish registers, we began the fitting window at the start of the outbreak window (see below); for LBoM data, we started immediately after the last observation before the peak for which the mortality rate was of the peak mortality rate. For two epidemics for which the data are particularly noisy (1593 and 1625), we were unable to find any simple heuristic that would automatically choose a good fitting window and resorted to choosing the fitting window by eye.

The fitting windows used for all our fits are listed in . Although we experimented with different methods for determining the fitting window (e.g., starting the window from the observation after the last local minimum before the peak or based on a threshold of of the peak mortality rate), we emphasize that our choices were entirely determined by seeking good fits (judged by visual inspection) and not by the pattern of estimated initial growth rates. Moreover, our conclusion that later plagues were much faster than the earlier plagues is robust to choices we made about window selection.

Outbreak Years and Outbreak Windows.

We refer to the calendar year in which an epidemic begins as the outbreak year. For each epidemic, we define an outbreak window that is the maximum time interval to be considered for fitting. For the 14th century epidemics, we defined the outbreak window to be the full outbreak year, except for outbreak year 1348 when the epidemic spanned both 1348 and 1349 (so we took the outbreak window to be the full length of both calendar years). For the late-epoch epidemics, we defined the outbreak windows as sequences of consecutive weeks during which at least one plague death was listed in the LBoM ().

Data and Code Availability.

Our epigrowthfit R package for estimating initial growth rates, which includes all of the data used in this paper, is available in Github at https://github.com/davidearn/epigrowthfit. The data can also be obtained from the International Infectious Disease Data Archive (http://iidda.mcmaster.ca). The additional code required to reproduce all of our results is available in Github at https://github.com/davidearn/plague_growth.