Archaeological Radiocarbon Database.

The publicly available dataset includes more than 13,658 archaeological radiocarbon dates comprising 2,759 sites for Mesolithic and Neolithic Europe, ca. 10–3.5 kya (65). Each site is assigned to a geographic region based on physiography, perceived spatial clustering, and prior knowledge of prehistoric cultural regions. Problematic dates including outliers and those with ambiguous stratigraphic associations are discarded. We select a subset of nine regions (Table S2) for the EWS analysis because they provide the clearest qualitative and quantitative shifts from high to low growth regimes–that is, population booms and busts (17), raising the possibility of EWSs leading up to the bust events.

To minimize the effects of radiocarbon sampling bias within and among sites, we use the following phase-based summing procedure (22⇓–24). Each date for a given site is assigned to a temporal site-phase and multiple phases are defined at a site when the minimum pairwise distances among all radiocarbon dates exceed 200¹⁴C y. One ¹⁴C assay is then randomly selected from each site phase and calibrated {1} using the Bchron R package (67) and IntCal13 atmospheric correction curve (68). The resultant probability distributions are summed to unity. This procedure is repeated 20,000 times and each SPD is averaged {2}. The resultant SPDs are binned in 10-y periods. As a final step, we correct each probability density curve for radiocarbon calibration effects (66) and taphonomic bias (21) by estimating a null model using a Monte Carlo routine and subtracting it from the observed curve to generate a corrected SPD curve {9}.

SPD Null Model Generation.

To explore the possibility of observing EWSs when they did not actually exist in the systemic context (type I error), we construct a null model that simulates monotonic population growth, archaeological sampling, atmospheric effects on radiocarbon proxies, and taphonomic decay of radiocarbon proxies (Fig. S3). We first fit a generalized linear model (GLM) to the empirical SPD {3}. This simultaneously models exponential growth of populations and decay of archaeological carbon (17). Second, to model sampling, we randomly draw a probabilistic sample from a uniform distribution using the fitted GLM for probability weights and a sample size equivalent to the empirical sample size {4}. Third, to model the effects of the radiocarbon calibration process, we “un-calibrate” each simulated date by passing it backward through the IntCal13 radiocarbon calibration curve {5}. We then calibrate the uncalibrated values using them as ¹⁴C means with SE values sampled without replacement from the SE values in the empirical sample {6} and sum the results {7}. This procedure is iterated 20,000 times and all calibrated probability distributions are summed to unity to produce a null SPD model {8}.

Calculating EWS Metrics Using SPD.

The methods for calculating EWS statistics generally follow those outlined in ref. 69. For each of the nine population trajectories, we isolate a subset of the time series from the beginning of the population boom to the end, or the collapse date {10}. The beginning date takes into account both qualitative evidence of growth as observed in the SPDs and the estimated date of earliest local archaeological evidence for agriculture. The date of collapse in each region is determined qualitatively by inspecting the SPD plots. The SPD curves may include edge effects due to decreased sampling in later periods, however the collapse dates are earlier in time so this does not affect the EWS analysis results. Next, we detrend the time series to remove linear structure in the time series by fitting a LOESS model and calculating residuals. The LOESS fitting procedure uses local regression to predict SPD values along the x axis (λ=2; α=0.15), and residuals are calculated as the difference between the predicted and observed SPD values. We then calculate autocorrelation [AR (1)], variance (σ), and skewness (1λ) values for the detrended data using a sliding time-window {12}. The first window includes the first 25% of points in the time series and then iteratively slides forward one decade at a time until reaching the defined collapse date. The size of this sliding window follows the recommendation of (69) and considers the scale of human demographic processes. Additionally, we conduct a sensitivity analysis (Dataset S1) for each region to evaluate the effect of sliding window size and the length of the time series on the significance of the EWS statistics (70). Autocorrelation within the sliding window is one of the simplest measurements of critical slowing down. AR (1) measures the degree of correlation between pairs of temporally adjacent measurements of state variable values (i.e., SPDt and SPDt+1). As the sliding window approaches the collapse point, increasing AR (1) values indicate the SPD values are becoming increasingly similar, which is an EWS for decreasing system resilience. This measurement is made here by fitting a lag-1 [AR (1)] ordinary least-squares regression model to each window of data and estimating the autoregression coefficient. More precisely, temporal autocorrelation is found as follows: zt+1=α1zt+εt, where α is an autoregressive coefficient and εt is Gaussian white noise. Variance, reported here as a SD (σ), quantifies the spread of values, which may indicate a system is approaching a regime shift. σ is found as follows: σ=[1/(n−1)]∑t=1n(zt−μ)2, where z is SPD density. Skewness (1γ) measures directional bias in the data as the standardized third moment around the mean of a distribution, or 1γ=[(1/n)∑t=1n(zt−μ)3]/(1/n)∑t=1n(zt−μ)2 (69). The shape of the distribution of system indicator variables can change in systems that are undergoing CSD because as the system slows down, it will hover around a critical transition point. These values can also include ones that indicate alternative system states, a process sometimes called flickering (32). This ambivalence in the system trajectory can yield increasingly skewed distributions with significant skewness statistics in either the positive or negative direction, depending on whether a system is moving toward larger or smaller values.

Once the summary statistics are calculated for each time increment, we qualitatively assess the data trends for EWSs (Table S1); if a given system exhibits EWSs of a regime shift (collapse), theory predicts increasing autocorrelation and variance and decreasing skewness. Next, we assess these predictions for each region statistically using the Kendall’s tau (τo) rank correlation test for the entire time series subset. The τ statistic indicates the direction of the EWS statistic as either increasing (0<τ<1) or decreasing (−1<τ<0). We also calculate an optimized tau value (τ[o]) by finding the maximum (or minimum) τ value from time series of varying durations that iteratively extend back from the collapse date beyond a minimum temporal window of 100 y before the collapse. This optimization routine is necessary for two reasons. First, we do not know a priori when resilience began to decline in these systems. The time series data under consideration may be much longer than the period of declining resilience. Second, the EWS:noise ratio is expected to be low and may preferentially obscure signals that are necessarily weaker earlier in time. Optimized tau values are calculated separately for each EWS statistic because we currently have no reason to expect different EWS metrics to change at proportionate rates and thus emerge through background noise at the same point in time. Rather, once a system begins to experience declining resilience, different EWS may change at different rates and thus emerge through background noise at different points in time. The EWS results for the Neolithic Paris Basin illustrate a situation in which (i) the quantitative method as applied to the full range of data are unable to capture qualitatively distinct AR (1) and σ EWS leading up to collapse, and (ii) AR (1) and σ trends indicate differential sensitivities to the empirical SPD. Conversely, the optimization routine identifies these trends, and, therefore, we report the τ[o] in Fig. 2 and Table 1. We follow standard approaches and report the probability values associated with τo {13} and τ[o] {14} using single-tailed tests, Ps and P[s], respectively {13, 14}. For completeness, we report τo in Table S3 and both τo and τ[o] in Dataset S1.

Recognizing the potential for spurious results given the complexity of these archaeological cases, we develop a rigorous statistical test to examine the effects of confounding factors. To do this, we analyze the SPD null model. We calculate EWS statistics for each of the 20,000 simulated density curves and combine them into distributions represented as {τe} and {τ[e]} {15}. Cumulative density functions (CDFs) are computed {16, 17}, and the probabilities of τo and τ[o] are determined for each EWS statistic, given a null model of exogenous processes {18, 19}. Thus, we define two “robust” probability values using single-sided tests, Pn and P[n] as F(x)=Pr(X≥x), where F(x) is the CDF for {τe} and {τ[e]}, respectively.

Last, we examine the probability of obtaining the observed ratio of statistically significant:insignificant results. We apply a Fisher exact test, which involves classifying the probability values for each EWS statistic of each region as either passing or failing given a one-tailed probability test Pr(X≤0.1) to determine the likelihood that the significant EWS statistics τo and τ[o] could have occurred by chance. Given the sample size of nine cases, we use a 0.1 significance level in this classification. The number of observed cases of nine that pass or fail a given EWS test is given using an expected pass:fail ratio of 1:8. We report Fisher exact test results for P[s] and P[n] in Table 2 and for Ps and Pn in Table S4.

A Regime-Shift Model for Human Settlement Densities at Large Spatial and Temporal Scales.

The model shown in Fig. 1 uses a simple step function with a high value step decreasing from 0.4% to 0.2% and a low value step at −0.4% growth. The upper limit of 0.4% growth reflects the upper end of expected annual population growth rates in early agricultural populations (12). To simulate increasing variance leading up to collapse, a randomization function is applied to the upper step values in which the mean is zero and the SD increases from 0.1% to 0.2%. The variance in the lower step is set to be constant at 0.2%. The growth rate model is then used to simulate an SPD curve by assuming a starting population of 500 individuals and iteratively increasing the population size by the modeled growth rate value at a given time step. We model the effects of the radiocarbon date calibration process, including atmospheric effects and isotopic counting error, by back-calibrating a uniform sample of dates 7–5.5 kya (n = 1,000), adding error sampled from a Poisson distribution (λ=50), and calibrating these simulated dates using the IntCal13 calibration curve (68). To model the effects of archaeological sampling, a random sample of 1,000 dates is drawn from the simulated population. The modified population curve is finally summed by the same binning procedure used on the empirical data. The result is a simulated SPD that results from a theoretical fold bifurcation in population growth rates, and we calculate the variance EWS metric on the SPD data assuming a collapse date at 6,225 cal. BP. We use a 77-y sliding window (∼10% of the 775-y duration; see ref. 69). Finally, we calculate τo and test its significance using Ps and Pn. We note that this model is developed to illustrate that EWSs are detectable in archaeological SPDs despite the potentially confounding factors discussed in the main text. The model explores a number of interactions, although it is not exhaustive. Other interactions, such as the relative strength of EWSs that would be required for detecting EWSs at varying points in the calibration curve or the possibility of CSD due to biological population growth may require additional investigation.

Wavelet Analysis.

A univariate continuous wavelet transform analysis (71) is run using the R biwavelet package (72). For each region, we analyze a corrected and detrended SPD and plot two 95% significance contour intervals (Dataset S1). The significance contour plotted in white uses the built-in time-averaged filters. Plotting the second significance contour involves performing univariate wavelet transforms and calculating the 95% quantile on the distribution of power matrices from the simulated SPDs from the null model. The power matrix from the observed SPD is then divided by the quantiles from the SPD null simulations, and the resulting 95% significance contour is plotted in black.

Early Warning Signals of Collapse in Neolithic Paris Basin.

Here we briefly step through a sample case to illustrate the procedure and interpretation of results. Fig. S1 illustrates how EWSs are detected in the Neolithic Paris Basin. Fig. S1B shows the complete SPD curve including the region of the time series isolated for EWS analysis. The SPD curve serves as a temporal proxy for population change throughout the Paris Basin Neolithic period. Following previous studies (17, 20, 26), we assign a regime-shift date, or “bust” point, at which time the system shifted from population growth to decline. For the Paris Basin example, the bust point was previously identified at 6,225 cal. BP. We then qualitatively define a point-of-beginning for the Neolithic boom. In most cases, this is simply the point at which agriculture is thought to have entered the region. However, if the SPD exhibits growth before the currently accepted agricultural start date, we assign an earlier point-of-beginning. This decision is justifiable given that archaeological estimates of agricultural origins reflect minimum possible dates due to sample limitations. For the Paris Basin, we assign the beginning of the population boom to 7,500 cal. BP, which precedes the currently accepted agricultural origins date by 300 y. Thus, the Paris Basin population boom isolated for EWS analysis ranges from 7,500 to 6,225 cal. BP.

Fig. S1A shows the extracted time series, the time-series trend model, and a schematic representation of the EWS sliding window, all of which are used in the derivation of EWS metrics. Fig. S1 C, E, and G shows the results of the EWS analysis for the three statistical indicators considered here, including lag-1 autocorrelation [AR (1)], SD (σ), and skewness (1γ). AR (1) and σ are expected to increase (+), whereas skewness is expected to decrease (−) in advance of regime shift (37). In the case of Neolithic Paris Basin, qualitative increases in AR (1) (+) and σ (+) are apparent after 6,800 and 6,300 cal. BP, respectively. In contrast, 1γ does not exhibit clear directionality, although some decrease may be present. Kendall’s tau (τ[o]) statistics are used to quantify directionality and to assess statistical significance. The τ[o] statistics support the qualitative observations at P[s]<0.01. That is, AR (1), σ, and 1γ exhibit significant departures from zero.

These statistics tentatively support the presence of CSD in the empirical data, but they cannot distinguish CSD due to systemic changes in human demography from CSD due to exogenous factors such as sampling error, atmospheric effects on radiocarbon, and taphonomic biases. Fig. S1 D, F, and H shows the cumulative density distributions [ECDF(τ[e])] for trend statistics generated from a null model simulating the effects of these confounding factors and the probability of obtaining the empirical values given the null model (Materials and Methods). The probability values for the τ[o] trend statistics for AR (1) and σ reveal significant departures from the null model (P[n]<0.1), whereas 1γ does not. Therefore, the CSD trends are more likely to reflect demographic processes rather than exogenous factors. We conclude that CSD preceded the Neolithic Paris Basin collapse that occurred around 6,225 cal. BP. In contrast, the skewness results are inconclusive.

Theoretical Considerations Related to SPD-Based EWSs Among Human Societies.

Although regime shifts and CSD may be expected, direct observation of such dynamics at requisite spatial and temporal scales is difficult. Archaeological data provide us with long-term, spatially extensive proxies for human demography (4). The SPD method in this study does not directly measure human population levels but rather approximates the number of inhabited settlements in a given spatiotemporal unit, which is taken to index human population levels after applying appropriate corrections. Thus, SPDs directly measure how settlement densities changed through time, but they do not directly measure population change. However, recent work has shown strong correlation between SPDs and a more direct demographic proxy based on skeletal measurements known as the juvenility index; this study also compared settlement plans and ethnographic censuses from contemporary farmers and foragers and found that differential settlement densities do not obscure demographic reconstruction using SPDs (26).

Because CSD is predicted to occur in human populations, a theoretical question that arises is how CSD should manifest in human settlement change, which may involve distinct dynamics. To explain how human settlement dynamics could come to reflect CSD in human populations, we hypothesize the relationships between archaeological settlement patterns in high and low growth regimes. In a high growth regime, such as occurred immediately after the introduction of agriculture in Neolithic Europe, there would have been a mix of young and old villages with new settlements forming from parent villages and rapid overall regional growth. In this regime, any abandoned settlements would have been rapidly repopulated because of high regional growth rates and high demand for resources. In a low growth regime, settlements would have been older on average with fewer new settlements appearing. Demographic variance and system shocks such as warfare, conflict, disease, or stochastic environmental effects would have had significant impacts, and failed settlements would have rarely repopulated because of the low overall demand for resources.