SINE QUA NON: INFERRING KODJADERMEN-GUMELNIȚA-KARANOVO VI POPULATION DYNAMICS FROM AGGREGATED PROBABILITY DISTRIBUTIONS OF RADIOCARBON DATES

Abstract

Past human population dynamics play a key role in integrated models of understanding socio-ecological change over time. However, little analysis on this issue has been carried out for the prehistoric societies in the Lower Danube and Eastern Balkans area. Here, we use summed probability distributions of radiocarbon dates to investigate potential regional and local variation population dynamics. Our study adopts a formal model-testing approach to the fifth millennium BC archaeological radiocarbon record, performing a region-wide, comparative analysis of the demographic trajectories of the area along lower Danube River. We follow the current framework of theoretical models of population growth and perform global and regional significance and spatial permutation tests on the data. Specifically, we investigate whether populations on both sides of the Danube follow a logistic pattern of steady growth, followed by a major decline over time. Finally, our analysis of local-scale growth investigates whether considerable heterogeneity or homogeneity within the region may be observed over the time span considered here. The results show both similarities and differences in the population trends across the area. Our findings are showcased in relation to the cultural characteristics of the region’s 5th millennium BC societies, and future research directions are also suggested.

INTRODUCTION

The so-called “Kodjadermen-Gumelniţa-Karanovo VI cultural complex” (KGK-VI) is a component of the Southeastern Eneolithic Block (SEB), defined based on material culture by the cultural-historical archaeologists from Bulgaria and Romania. It was identified within a broad region of the Eastern Balkans and Lower Danube Valley, delimited to the north by the Carpathians mountains, but also the steps area from north of the Danube Delta to the east by the Black Sea and reaching as far as the Rhodope mountains and Olt River to the west and the Aegean in the south. The general chronological position of KGK-VI evolution is placed during most of the 5th millennium BC. It is characterized by the widespread development of tell settlements, the emergence of extramural cemeteries (in some cases including wealthy graves), the advent and development of copper and gold metallurgy, along with consistent changes in lithic and ceramic technologies (e.g., graphite and gold pottery painting, exploitation, and processing of various pigments, etc.). Local variants emerged from a common, prior cultural background (e.g., “Varna Culture” on the western coast of the Black Sea) (Todorova and Zhelyaskova 1978; Todorova 1986; Marinescu-Bîlcu 2001; Popovici 2010; Petrova 2016; Lazăr et al. 2018; Chapman 2020) and also in areas of interaction with other cultures (e.g., known as Stoicani-Aldeni/Aldeni II cultural group in north-eastern Wallachia) at the intersection between KGK-VI and Cucuteni-Tryplie communities (Frînculeasa 2016), or in previously unoccupied territories (e.g., known as Bolgrad/Bolgrad-Aldeni group) in the region north of the Danube Delta (Todorova and Zhelyaskova 1978; Todorova 1986; Frînculeasa 2016), which marks the most northern and eastern extension of the KGK-VI.

Moreover, the increased mobility of the KGK-VI population is highlighted by a vast exchange network of raw materials and goods. Thus, pottery painted with black, white, and red pigments, graphite or even with gold, along with metal or flint items (e.g., super-blades), adornments made of Mediterranean shells (e.g., Spondylus sp., Dentalium sp., etc.), along with non-local minerals and pigments (e.g., malachite, marble, carnelian, agate, hematite, kaolin, etc.) are common discoveries in the KGK-VI sites (e.g., tells, flat settlements, off-tells, and cemeteries) located in areas that usually lack all of these raw materials (Bailey 2000; Todorova 2003; Anthony et al. 2010; Popovici 2010). All of this reflects the adoption of new lifeways, economic development and a new way of environmental exploitation and control, in parallel with the development of complex and stratified society and the emergence of elites, as demonstrated by the wealthy graves from Varna I cemetery or other domestic/funerary discoveries (Bailey 2000; Chapman et al. 2007; Anthony et al. 2010; Slavchev 2010; Chapman 2020). Consequently, starting from that flourishing development of the SEB human groups, without any parallel in the past, some scholars described the period as the “Ex Balcanae Lux” phenomenon (Todorova and Zhelyaskova 1978; Sterud et al. 1984), the “Climax Copper Age” (Chapman et al. 2007) or the “Golden 5th Millennium BC” (Boyadzhiev and Terzijska-Ignatova 2011), in order to highlight the extraordinary progress of the Neolithic communities.

Within this large region, different human groups have developed specific material culture signatures, or archaeological “cultures,” during this period. Recent aDNA investigations in southeastern Europe have demonstrated that populations responsible for the above mentioned “cultures” of the second half of the 6th millennium and 5th millennium BC, KGK-VI included, have similar genetic features and common ancestry due to their common origin in southwestern Anatolia (Hervella et al. 2015). Along with the Anatolian Neolithic ancestry, those populations showcase sporadic evidence of steppe-related ancestry (specimens in Varna I and Smyadovo cemeteries in Bulgaria), but also a consistent hunter-gatherers related ancestry (some resilient native Mesolithic groups from the target area), which indicates a complex population structure and admixture, with several genetic components (Mathieson et al. 2018).

The last decade has seen major advances in developing theoretical, analytical, and methodological instruments, concerning the understanding of demographic change out of large datasets of archaeological radiocarbon (¹⁴C) dates, in different parts of the world and encompassing large temporal spans of the Late Pleistocene and Holocene (Shennan et al. 2013; Downey et al. 2014; Timpson et al. 2014; Crema et al. 2016; Downey et al. 2016; Bevan et al. 2017; Barton et al. 2018; Riris 2018; Roberts et al. 2018; Timpson et al. 2020; Crema and Shoda 2021). Some of these studies have also been focused on understanding the population dispersals and change during the Neolithic of Southeastern Europe, although mostly on the early Neolithic of the region (Porčić and Nikolić 2016; Blagojević et al. 2017; Harper 2019; Vrhovnik 2019; Porčić et al. 2021; Vander Linden and Silva 2021).

The current study aims at contributing to the understanding of the population geo-temporal dynamics of one of the archaeological signals of the Chalcolithic in Eastern Europe, namely KGK-VI, which reflects the maximum point of development of human communities that lived a Neolithic way of life in this part of Europe. More specifically, our research explores (1) the nature of population trajectories within the chosen region, on both sides of the Danube and (2) the extent to which KGK-VI differs north and south of Danube. The analyses used in the study allow for local and global tests of significance to be performed and regional population histories to be compared through the comparison of empirical and simulated summed probability distributions (SPDs) of radiocarbon dates (see Bevan et al. 2017; Crema and Bevan 2021 for details).

DATA

In the current study 440 radiocarbon dates from both Romania and Bulgaria ascribed to the KGK-VI were used (Figure 1). The data were acquired from publications, gray literature, unpublished sources, and from on-line databases (Reingruber and Thissen 2017; Weninger et al. 2020) (e.g., http://www.14sea.org/index.html and https://www.academia.edu/40774947/CalPal_Holocene_Palaeolithic_14C_Database). The available radiocarbon dates from the above-mentioned sources were compiled in a database (n = 257 from Romania and n = 183 from Bulgaria), grouped into 135 bins for analysis (see explanation below of binning), recovered from 59 sites (n = 32 in Romania and n = 27 in Bulgaria). The database includes contextual information such as site name, site id, site recovery context, “culture” and phase (where available), region, country, laboratory number, the uncalibrated date and uncertainty, and geographical coordinates (longitude and latitude) (Supplementary Table 1).

Figure 1 Map of radiocarbon distribution data set. Legend: 1-Akladi Cheiri; 2-Bikovo; 3-Čardako-Slatino; 4-Djakovo; 5-Dolnoslav; 6-Drama-Merdžumekja; 7-Durankulak; 8-Ezero; 9-Goljamo Delčevo; 10-Hotnica; 11-Junacite; 12-Karnobat; 13-Košarna; 14-Omurtag; 15-Orlitsa; 16-Ovčarovo; 17-Povelyanovo; 18-Smjadovo; 19-Sušina; 20-Tatul; 21-Tell Azmak; 22-Tell Karanovo; 23-Tell Russe; 24-Varhari; 25-Varna1; 26-Varna2; 27-Varna3; 28-Baia Boruz Tell; 29-Popina Blagodeasca; 30-Bordușani; 31-Carcaliu; 32-CăscioareleOstrovel; 33-Cunești; 34-Dambul lui Haralambie; 35-Gumelnița-terrasse; 36-Gumelnița-tell; 37-Hârșova; 38-Lișcoteanca-Movila Olarului; 39-Lunca; 40-Luncavița; 41-Mălăieștii de Jos; 42-Măriuța-C; 43-Măriuța-T; 44-Navodari; 45-Niculițel; 46-Orbeasca Sus; 47-Panduru; 48-Pietrele; 49-Seciu; 50-Șeinoiu; 51-SultanaGhețărie; 52-Sultana-Malu Roșu-terrasse; 53-Sultana-Malu Roșu-tell; 54-Taraschina; 55-Taraschina_2; 56-Urlați; 57-Vărăști; 58-Vitănești; 59-Vlădiceasca.

The dates selected to use in the current investigation were obtained from samples from various materials including, wood, charcoal, seeds, and bones (herbivores and humans). The dates based on shell, fish and carnivore samples have not been included to avoid potential reservoir effects issues. We also applied a cleaning protocol and excluded all dates with large uncertainties in order to remove any potential spurious dates. However, most of our dates have very low uncertainties, with a mean of 43 and a median of 37 years. Recent studies, dedicated to similar research agenda, have also shown that an overly strict data cleaning and the exclusive use of dates with very low uncertainties might potentially be just as damaging for this kind of analysis just as an uncritical data collection (Timpson et al. 2014, 2015; Vander Linden and Silva 2021).

METHODS

Our analysis was carried out with the R environment for statistical analysis (R Core Team 2020), and the rcarbon R-package for date calibration and SPD modeling (Bevan et al. 2020; Crema and Bevan 2021), using the Northern Hemisphere calibration curve (IntCal20) (Reimer et al. 2020), along with gstat (Pebesma and Graeler 2021), sp (Pebesma et al. 2022), and other R packages mentioned in the rmarkdown script that accompanies the study. Dataset, supplemental figures and R Markdown scripts needed for reproducing the results of our analysis are available at: https://zenodo.org/record/7587242#.Y9hI9S8Ro4c. The methods that we used in our analysis are based on previously developed quantitative analysis of SPDs by Shennan et al. (2013), further refined in several other recent studies (Downey et al. 2016; Brown and Crema 2019; Crema et al. 2016; Riris 2018; Timpson et al. 2014; Timpson et al. 2015), the non-parametric extension devised by Crema et al. (2016; see also Crema and Kobayashi 2020), and the spatial permutation test (Crema et al. 2017; Crema and Bevan 2021). The reader should refer to these bibliographical resources (and references therein) for more details on the methods and the concepts behind them.

Spatial Dispersal Analysis

We first analyzed the spatial structure of the dispersal of our ¹⁴C distribution (Figure 2a,b). We considered this as being the first step in our analysis for a better understanding of the spatial dynamics of the dispersal of the KGK-VI to assess its correlation (potentially) with demography. Based on Hengl (2006) equation for establishing the right pixel size for a grid, and our regional context we superimposed a 23 × 23 km grid on the area, with each grid cell covering approximately 460 km². The analysis was done in R statistical environment with the gstat (Pebesma and Graeler 2021) and sp (Pebesma et al. 2022) packages. We divided the study area into grid cells and hexagon cell shapes were chosen, given their shape being the closest to a circle and the easy to use in a tessellation. Their minimal edge effects as well as the identical neighboring cells and having the same distance between centers for all the neighbors, make them particularly suitable for our analysis (see also Vrhovnik 2019).

Figure 2 (a) Number of radiocarbon dates per grid cell. Values are log10 scaled. (b) The earliest appearance of the KGK-VI settlements in the grid cells, consisting of grid cell centroids with the date for the beginning of the KGK-VI occupation. Gridded area (white hexagons) represents the KGK-VI area with dated sites.

Thus, the study area is covered with 477 grid cells, of which only 47 grid cells are occupied with sites, forming several clusters, and a patchy distribution of samples. The number of radiocarbon dates per grid cell varies from one (seen in 11 grid cells) to a maximum of 67 (seen in one cell) with a median value of 3 dates per grid cell and third quartile at 10.5 dates per grid cell. The distribution of dates per cell is shown in Figure 2a.

For each grid cell, we then calculated a normalized summed calibrated radiocarbon probability distribution. To calculate the calendar age ranges, highest probability density was used, and these are the shortest ranges that include 95% of the probability in the summed probability density function. As such, the starting date of the KGK-VI in a particular grid cell was taken to be the lower 95% range endpoint date. These estimated starting dates are shown in Figure 2b. Consequently, these dates were to estimate the spread of the KGK-VI across the area. Grid cells with only one radiocarbon date were excluded from the interpolation.

Dates Binning, Calibration, and SPDs Production

To mitigate potential issues due to the differences in intensity of sampling, radiocarbon dates are combined in 100-year bins within each archaeological context (e.g., horizontal and vertical provenience units) so that the intensively sampled sites/areas are not overrepresented and cause artificial spikes in the observed SPDs. When multiple dates are present from a single site, they are aggregated within each archaeological context and when their distance in ¹⁴C years is less than 100 years. The procedure applies a hierarchical cluster analysis using the complete linkage method and a cut-off value of 100 years to separate the observations. Although our selection of 100 years for binning is arbitrary, we performed a bin sensitivity analysis (Supplementary Figure 1), which shows this choice has no negative impact on the accuracy of results (all bin sizes fit within the 95% confidence simulated envelope—gray area) and is also well above the median error of 37 years in the dataset (our protocol follows those already applied in the literature; for more details (Bevan et al. 2017; Riris 2018; Crema and Bevan 2021).

Dates were calibrated using the Northern Hemisphere Radiocarbon Age Calibration Curve (IntCal20) (Reimer et al. 2020) and the rcarbon package (Bevan et al. 2020). Multiple dates within a bin are calibrated and summed “inside” the bin and subsequently divided by the number of dates so that each archaeological context contributes a single date distribution to the overall SPD. The probability distributions of the calibrated dates were summed over the entire KGK-VI period to produce empirically based SPDs using the entire data set, as well as for subsets for the two regions north and south of the Danube. Following detailed discussions in recent works (Weninger et al. 2015; see also details in Bevan et al. 2017), we have not normalized the post-calibration distribution of each date (that ensures it sums to 1 under the curve) before summation of multiple dates. This ensures the reduction of creation of abrupt spikes in the final summed probability distributions, there where the calibration curve is steep (Weninger et al. 2015; Bevan et al. 2017; Crema and Bevan 2021).

Model Testing

In order to differentiate SPD fluctuations that represent meaningful demographic change from those due to sampling error noise, we compare the observed SPDs against null models of simulated dates derived from hypothesized calendar age distributions of dates. These null (hypothesized) models assume increasing survival of datable radiocarbon material through time according to either an exponential or a logistic population growth. Only portions of the SPD curves that fall outside the 95% confidence interval (CI) of the null models are considered sufficiently significant demographic changes to be considered in our subsequent discussion. Model testing procedures are described in detail below.

We first evaluate the goodness-of-fit of the entire data set SPD, by first fitting the calibrated data to a generalized exponential model, with the help of modelTest function (Figure 3) (Timpson et al. 2014; Crema and Bevan 2021) (Supplemental Material for analysis in R). An exponential model had become common practice, as it is assumed to account for population growth with unlimited resources and taphonomic processes. We assessed whether the SPDs of the ¹⁴C dates for the entire region showed statistically relevant deviations when compared against the exponential model, following the procedure described in several other studies (Timpson et al. 2014; Bevan et al. 2017; Crema and Bevan 2021; Crema et al. 2016; Riris 2018; Palmisano et al. 2017; Roberts et al. 2018 and references therein).

Figure 3 Results of fitting and comparing the entire regional empirical SPD against the exponential null model of population growth. Monte Carlo 95% confidence null model gray envelope is based on 5000 runs. Observed SPD is shown with solid red line, while the positive and negative deviations from the null are marked in red and blue. (Please see online version for color figures.)

Null models were simulated for the entire KGK-VI region using assumptions of the exponential and logistic growth patterns in the following way, repeated for each type of model.

1. 1. An exponential/logistic model was generated for the entire time period and fit to the empirical SPD produced by the dates we compiled.

2. 2. A set of dates (equivalent in number to the bins used for the empirical SPD) was generated from the model and errors assigned randomly (within the range of empirical date errors).

3. 3. An SPD was generated from the model dates and errors.

4. 4. Steps 2–3 were repeated 5000 times to estimate the 95% CI around the model.

The empirical SPDs are then compared with the 95% CI around each of the two models (exponential and logistic). Portions of the observed regional SPD that fall outside the envelope were considered statistically significant local deviations above and below the null model (red and blue areas respectively). Following methods outlined by Timpson et al. (2014) a global p-value can then be calculated from the total area of the empirical SPD curve that falls outside the 95% CI of each null model.

To evaluate how well the exponential model fit to our empirical data we employed Akaike Information Criterion (AIC) to select the most parsimonious fitted model for the entire region (Sakamoto et al. 1986). AIC suggested that the use of different model might be a better fit for our context. As such, we decided to use the logistic growth model as the null model for comparison and discussion of results (Figure 4), following the same procedures outlined above for the exponential growth model. We also used an extension of the global p-value to evaluate the point-to-point differences along the SPD curve. This test compares the empirically observed difference to the distribution of differences in the SPD curve between two points in time against the distribution of expected values under the null hypothesis (Edinborough et al. 2017). We also deployed a more recent alternative to SPDs, Composite Kernel Density Estimates (CKDEs) (Brown 2017; McLaughlin 2019), which has the advantage of minimizing calibration artificial spikes, as well as, providing estimates of sampling and calibration-derived uncertainty over time (Figure 5). Based on the qualitative inspection of the SPDs, the CKDE curve and formal AIC test, we decided to also fit a suite of four composite models (see Rmarkdown document in the Zenodo repository) to the summed calibrated probability distribution for KGK-VI, representing potential demographic models (see Goldberg et al. 2016; Arroyo-Kalin and Riris 2021; de Souza and Riris 2021). Assessing for the goodness of fit of these models was achieved following the same procedure as outlined above.

Figure 4 Results of fitting and comparing the entire regional empirical SPD against the logistic null model of population growth. Monte Carlo 95% confidence null model gray envelope is based on 5000 runs. Observed SPD is shown with solid red line, while the positive and negative deviations from the null are marked in red and blue.

Figure 5 Bootstrapped composite kernel density estimate, suggesting a composite model with the breakpoint at approximately 4400 BC should be tested.

Permutation Tests—Regional Comparison

We are obviously interested in empirically testing the variation between the northern and southern regions of the KGK-VI (Danube River is used as a geographical divide). In order to achieve this, each region was compared to a null model assuming no spatial differences. This null model can be obtained by pooling the radiocarbon dates from across the entire region and simulating from the pooled SPD (Figure 6).

Figure 6 The KGK-VI empirical SPD record fitted to logistic (5000 BC–4400 BC) and exponential (4400 BC–5750 BC) models, with significance envelope derived from 5000 Monte Carlo simulations.

Crema et al. (2016) developed a permutation-based test to statistically compare two or more SPDs (Figure 7). The null hypothesis is generated by simulating multiple (e.g., 5000 here) SPDs whose dates are drawn randomly from both subregions (north and south of the Danube). These simulated SPDs are again combined, as above, to produce a 95% CI envelope against which the SPD from each region can be compared (see Crema et al. 2016; Crema and Bevan 2021).

Figure 7 Permutation test showing variation between regional population growth. Observed SPDs for each region are shown with a solid black line, while the dashed line represents the observed pan-regional SPD. Gray areas represent the 95% confidence envelope for the null model, red and blue bands represent areas where the observed SPD significantly positively (red) and negatively (blue) deviates from the pan regional null model.

This is a robust approach to inter-regional differences in the research intensity because the comparison is based on the shape of the SPD (the relative change in summed probabilities within each region) and not on differences in their absolute magnitudes. Maintaining the observed number of bins for each region and comparing population trajectories rather than absolute differences in density, the permutation test bypasses the problem. Moreover, sample size is taken into account in the width of the 95% CI envelope. Significant negative (or positive) deviations of the SPD in one region does not necessarily imply a lower (or higher) absolute population density, but that the drop in the proxy within the dynamics of that region was significantly stronger compared to rest of the data.

Spatial Permutation Test

The spatial permutation test is an extension of the permutation test described above, having the virtue of allowing for the assessment of variation without the imposition of a priori regions of analysis. The steps involved in the spatial permutation test are described in detail by Crema et al. (2017; see also Crema and Bevan 2021). Below are summarized the steps involved in the spatial permutation analytical protocol.

1. 1. Produce local SPDs for each site with dates combining the date distribution at the site with date distributions at neighboring sites weighted as a function of their distance from the site. We selected a neighborhood radius of 100 km following a sensitivity analysis of different radii (see Supplemental Figures 2–3).

2. 2. Divide the KGK-VI temporal span of 5050–3800 BC into equal transition blocks (e.g., time slices), 250 years each, which, given the length of our time span, we consider relevant in our endeavor to detect long term regional growth patterns (Figure 8).

   

   Figure 8 Observed rate of growth at each transition computed from the SPD. I: 5050–4800 to 4800–4550; II: 4800–4550 to 4550–4300; III: 4550–4300 to 4300–4050; IV: 4300–4050 to 4050–3800 BC.

3. 3. Calculate the overall growth rate (change in the SPD curve) between each temporal transition block and the subsequent one.

This allowed us to evaluate spatial patterns of demographic growth and decline in two ways. We compared the growth rates calculated for each local SPD with the overall growth rate for each transition block (Figure 9). For each transition block, we also compared the growth rates of local SPDs at each site with a simulated model generated by repeatedly (10,000 iterations) randomly shuffling the local SPDs spatially across all site locations and combining the growth rates of the shuffled SPDs at each site. This allowed us to identify “hot” and “cold” spots (areas of significance), defined as areas where the local growth exceeds the growth observed in the simulation (Figure 10). Following methods discussed in Crema et al. (2017), two measures of significance are produced in the course of the spatial permutation test. p-values are measures of significance between observed local growth and simulated growth rates. However, the use of multiple testing approach, increases the potential for compounding false positive results (e.g., some local SPDs will be higher or lower than the theoretical expectation by chance alone). A more robust q-values test is therefore also computed by adjusting p-values to account for false positive discovery rate. Thus, a p-value of 0.05, implies that 5% of the tests will result in false positives, a q-value of 0.05 means that 5% of the results that have a q-values less than 0.05 are false positives (see Crema et al. 2017 for further details).

Figure 9 Local geometric growth rate for each transition block. I: 5050–4800 to 4800–4550; II: 4800–4550 to 4550–4300; III: 4550–4300 to 4300–4050; IV: 4300–4050 to 4050–3800 BC; V: Geographical reference map for the local geometric growth rate analysis, shown in transitional blocks I-IV.

Figure 10 Spatial permutation test showing where growth is significantly higher or lower than the null for each transition block. I:5050–4800 to 4800–4550; II:4800–4550 to 4550–4300; III:4550–4300 to 4300–4050; IV:4300–4050 to 4050–3800 BC; V: Geographical reference map for the spatial permutation test analysis, shown in transitional blocks I-IV. Significance is shown in terms of q-values (more robust against false positives) and p-values.

RESULTS

Spatial Dispersal Analysis

Figure 2(a–b) displays the results of spatial interpolation based on the analysis of the earliest appearance of the KGK-VI settlements in the grid cells. Early presence of KGK-VI across the region, between ∼4900 BC–4700 BC, although rather patchy, seems to be documented in several spots across the region, but mostly in the central and southern KGK-VI area, along major river valleys and their surroundings (see also in conjunction with Figure 1). The interpolation results thus suggest a fast emergence of this archaeological signature on both sides of Danube River, in its main core areas, in central-northeastern Bulgaria and southern Romania, extending in southeastern Romania, southern and southwestern Bulgaria (Figure 2b). There also seems to be a lag in KGK dispersal, roughly around 4600–4550 BC, or a weaker signal, that could be also due to a smaller sample of dated sites. From these core area KGK-VI dispersed in most of the region by 4400 BC. This analysis also suggests that potentially two more clusters, seems to emerge, one in south-southwestern Bulgaria and the other in southeastern Romania. Given the number of sites with radiocarbon dates, available for our study, we have not split the database into more clusters, as this might be the result of the current sample size and might have created artificial patterns in the SPD analysis and its model testing.

Model Testing

Figure 3 presents three noticeable fluctuations in the empirical SPD, from the exponential null model. The first of those is a significant negative departure from the exponential null model between 4895 BC–4661 BC (marked in blue). This is followed by a long significant positive deviation from the null between 4560 BC and 4281 BC (marked in red), with its peak at ∼4408 BC. A point-to-point statistical test for changes in the shape of the SPD, suggests that the difference between the peak in the empirical SPD and a first major trough at 4175 BC are statistically significant at a p = 0.0004 level. After this highest point the SPD drops rapidly, while remaining above the 95% CI null envelope, until ∼4300 BC. Subsequently, the SPD remains within the null model 95% CI envelope until dropping below this limit after 4200 BC until the end of the sequence, thus marking the second statistically significant negative deviation from the null model. The data, overall, exhibit significant deviations from the exponential model at a global level of p = 0.0002.

The result provided by the AIC test, comparing the fitted exponential and logistic null models indicated the logistic null model as more parsimonious (ΔAIC = 235, 42).

We reran the modelTest function with the logistic null model, and its result are shown in Figure 4. The SPD deviates statistically from the logistic model as well, with an overall p = 0.0002. A short positive deviation from the null model at the beginning of the SPD is most likely a statistical artifact of the method, which fits the model at a later start date. Inspecting Figure 4 we observe a better agreement between the fitted model (dashed line) and the 95% CI envelope, although toward the end, the envelope departs from the fitted line, and is more in line with the pattern of the empirical SPD. At the same time, the empirical SPD is in good agreement with both the null model envelope and the fitted model for about a half of the time range. The beginning and early stages of KGK-VI follow a steady growth pattern within the limits of the null model, until a rapid rise that starts approximately around 4600 BC and significantly deviates positively from the null and fitted line model, from 4553 BC until 4303 BC, reaching the peak around 4400 BC. The result of the point-to-point test that we have run between the highest peak in the SPD and the first major subsequent trough in the SPD around 4180 BC, also proved to be statistically significant at p = 0.0004 level. A final major deviation from the logistic model is represented by a long significant departure from the null, marking the final stages of the KGK-VI.

The bootstrapped CKDE analysis (Figure 5), as well as the shape of the empirical SPDs, and the significantly departure from both the exponential and logistic growth models, confirmed the inadequacy of both models to explain the KGK-VI data. We found that the best fit model (a Logistic-Exponential model (Figure 6), adequately describes the observed SPD by combining a phase of logistic growth (∼5000–4400 cal BC), followed by an almost symmetrical, this time exponential decrease (∼4400–3800/3750 cal BC), with a global p-value of 0.5257 (see also R code in Supplemental Material for the analytical protocol). It is worth mentioning that the empirical SPD, highly matches the 95% confidence envelope of the composite null model and that, while neither local nor global deviation from the null emerges, the decreasing trend after 4400/4350 BC is continuous with no signs of positive peaking, until the end of the KGK-VI archaeological signature.

Permutation Test—Regional Comparison

The permutation approach compares the regional trajectories from the two major areas of the KGK-VI, namely north and south of Danube, with the model representing the general growth trends of the targeted region, as opposed to an idealized model. Figure 7 shows that both regions significantly deviate from the whole region null model at p = 0.001 for the southern area, and p = 0.001 for the northern area, respectively. The overall results, thus suggests regional differences relative to the timing and nature of population change.

The south of the Danube area (Figure 7 bottom) starts positively beyond the limits of the regional null model, and the pan-regional SPD (dashed line) as well, until around 4580 BC, growing faster positively deviating from the null. The growth continues afterward but within the limits of the null model and reaches its peak at ∼4500 BC, but faster than the pan-regional empirical SPD (e.g., south of the Danube demography follows a faster pattern than the entire region as a whole at this time). Demographic stability is maintained at the highest density between approximately 4500–4400 BC, followed by an abrupt declining until 4200 BC, and continuing at a steadier pace until the end of the sequence. Two significant negative deviations from the null occur in the southern SPD between 4359–4320 BC and 4243–4222 BC. The north of the Danube area (Figure 7 top), on the other hand, displays a quite different trend initially displaying a significant lag against the null model, until approximately 4550 BC, followed by a faster growth, but closely matching both the regional null model (gray shaded area), and the pan-regional empirical SPD. Unlike the southern area, the northern density peak is reached at the same time as the region wide SPD (dashed line). The maximum density is maintained for about 100 years, like in the southern area. The declining pattern in the northern SPD follows a similar path like the southern one. The difference here is that there are two significant positive deviations from the null between 4358–4320 BC and 4256–4204 BC, respectively, contemporary with the negative ones in the southern region. Most of the local SPD declining pattern, however, falls within the 95% CI of the pan-regional growth model.

Spatial Permutation Test

The final step in our analysis evaluates the spatial variation in population growth without an a priori regional classification but based on the geographical coordinates (see details in previous section). Figure 8 shows the pattern of the observed pan-regional geometric growth rate and in the four intervals across the five 250-year transition blocks at which growth was measured (I: 5050–4800 BC to 4800–4550 BC; II: 4801–4550 BC to 4550–4300 BC; III: 4550–4300 BC to 4300–4050 BC; IV: 4300–4050 BC to 4050–3800 BC). This shows that the pan-regional trend had an initial rapid growth followed by a declining reduction of growth that extends even below zero once the carrying capacity is reached.

Evaluation of the local growth rates across the four temporal blocks (see Figure 9), shows some regional variation in actual growth rates, especially between the transitional blocks I-II and III-IV. To test the significance of actual growth rates and overcome potential issues caused by the calibration, we compared the local growth patterns to a null model based on regional wide growth trends (Figure 10). This allowed for the identification of significant positive or negative deviations from the overall trend. For the first interval of 5050–4800 BC to 4800–4550 BC (Figure 9), all KGK-VI sites from both sides of the Danube River present positive growth rates. However, as shown in Figure 10, most of the sites do not significantly positively deviate from the overall demographic trend, and only a few of them displaying negative deviation at p < 0.05. In the second temporal interval of 4800–4550 BC to 4500–4300 BC (Figure 9), while mostly displaying positive growth, some regional differences between north and south of Danube are apparent. Growth rates in the south range from 0.0004–0.004% annual rate, as opposed to 0.004–0.009% annual rate in northern area. The significance test (Figure 10) also shows obvious differences between the northern and southern KGK-VI areas. As such, the actual positive growth rates significantly depart from the regional trend in most of the northern region (q < 0.05), while the south trends in the opposite direction, with most of the area showing negative deviation from the overall pattern model, p < 0.05 level (see more details in e.g., Crema et al. 2017).

The third and fourth temporal intervals, 4550–4300 BC to 4300–4050, and 4300–4050 BC to 4050–3800 BC, respectively, show similar negative rates throughout the KGK-VI distribution area (Figure 9). However, as shown in Figure 10, these rates neither positively nor negatively deviate from the expectation of the pan-regional model, across the entire area.

DISCUSSION

We have developed in this research a set of radiocarbon-based demographic models of a well-known Chalcolithic archaeological signature in the Eastern Balkans and Lower Danube area, namely KGK-VI. This work contributes to other recent studies assessing Neolithic demographic dynamics in Southeastern Europe (Porčić et al. 2016; Blagojević et al. 2017; Vrhovnik 2019; Vander Linden and Silva 2021).

Using the spatial interpolation of the radiocarbon record we were able to show the spatiotemporal structure of the spread of the KGK-VI archaeological signature. It was marked by an early patchy manifestation in the main core area, as well as a few other regions both north and south of Danube. Thus, the analysis suggests a faster emergence of the KGK-VI archaeological signal in central-eastern, north-eastern areas of current Bulgaria, and spread across the entire studied area Although there is a slight lag in the expansion between south and north, the KGK-VI extends its entire occupation approximately simultaneously by about 4500 BC. However, given it rather rapid emergence and spread in some spots of the KGK-VI area, it is difficult, at this time, and it is beyond of this study aims, to determine whether this was due to a population dispersal, or to a rather demic-cultural combined model. This is something to be further modeled in other studies dedicated to this subject.

The positive rejection of the exponential model between 4560 BC and 4281 BC, demonstrates a higher-than-expected growth during this segment. The negative deviation from a null model in the early phase agrees with such a logistic growth trend. It may also be influenced by the slower pace in the northern area. The declining trend is also at a steeper and faster rate than expected, possibly indicating that population overshot regional carrying capacity. This agrees with demographic trends noted for the Neolithic observed in other areas of Europe (see below). The results of the regional permutation test display the regional differences between the north and south of Danube densities, however mostly in the details of their respective SPDs. As such, the northern area does show a lag in the local SPD, compared to the pan-regional null, while the southern area, displays a faster than expected pace toward reaching the peak in its density. However, looking at the shape of the curve of each region one can distinguish that their overall shape is quite similar and that their peak in summed densities fall within the 95% confidence pan-regional envelope, differing only in their position above (in the south) and partially below (in the north) the pan-regional SPD. The fact that the two regions show reciprocal patterns of demographic expansion and decline, suggest that we are potentially seeing some population shifting from the north to the south side of the Danube ca. 4600/4550 BC, and then some movement back the other way during the period of demographic decline after 4400/4350 BC.

Several studies of Early and Middle Holocene societies in other European regions, have shown “booms” followed by important “busts” (drops) in the SPDs (Shennan et al. 2013; Bernabeu Aubán et al. 2016; Downey et al. 2016; Bevan et al. 2017; Roberts et al. 2018; Palmisano et al. 2019), which are proposed to be related to long-term decline of Neolithic farming society. A similar pattern is shown in our study, when after a statistically significant “boom” and an equilibrium (of about 100 years—ca. 4450–4350 BC), a long period of higher than expected decline follows (about 350 years—ca. 4150–3800 BC), statistically significant against both exponential and logistic null models. However, we need a larger database to analyze SPDs for farming and climate and/or environmental co-variates, to help evaluate such potential explanations. Other potential explanations could be related to non-demographic aspects, such as changes in settlement patterns and in settlement size. Small more dispersed villages, during the early stages of KGK-VI, as well as a change from the larger complex tells, off-tell habitations, and flat settlement system from its pinnacle, toward much smaller more dispersed villages, may have increased the number of residential features that were dated, thus inflating the SPDs, and giving the appearance of higher population density, although it did not actually change. This is a research direction that deserves further investigation (see Porčić et al. 2021; Supplemental Material S1), as well as using complementary data, such as mortality juvenility index, and its correlation with the shape of the SPDs (see details in Downey et al. 2014).

CONCLUSION

In this exploratory study, we provided a possible reconstruction of the population dynamics of one of the well-known archaeological signals of the Chalcolithic in the Balkans and Lower Danube area, namely KGK-VI.

The current approach highlights the following conclusions:

1. 1. In the targeted geographical area, the human population density steadily grew for the first approximately 300 years following a logistic growth model until it reached the carrying capacity.

2. 2. The current study clearly shows the existence of at least two initial cores (north and south of the Danube) where what we call KGK-VI originally appeared in the early 5th millennium BC (4901–4751 cal BC) (Figure 2).

3. 3. The initial dispersion of the KGK-VI archaeological signal from each of the two cores, headed, for the first ca. 300 years, toward south/southwest on either side of the Danube, to disperse to the rest of the region by 4500–4400 BC, reaching the maximum development, documented until now for KGK-VI.

4. 4. The growth episode was followed by an equally fast start of the decline, beginning at ca. 4350 BC, without a significant recovery, until it vanished from the archaeological record toward the end of the 5th millennium BC–early 4th millennium BC.

5. 5. This decline was a continuous process, and it covers approx. 550 years (ca. 4350–3800 BC). The underlying causes behind both characteristics, especially relative to the declining trend, require more in-depth research. Therefore, the challenge for the future is to explore whether the “boom” and “bust” is related to deep-time historical contingencies under multiple factors.

6. 6. Our approach demonstrated the displacement vectors of KGK VI populations in the approached geographical area in correlation with the distinct temporal moments that marked them.

7. 7. Last, but not least, the present study demonstrates, once again (if necessary!), the fragility of the notion of archaeological “culture” (or cultural complex, group, aspect or facies) concept (-s), and especially that this Kossinian notion creates a major confusion between the notions of population, civilization and ceramic styles “fashion” (to which unfortunately we are still tributary). The entanglement created by these artificial notions in Balkan archaeology (and not only) led to excessive fragmentation of the prehistoric archaeological landscape. It prevented proper understanding (sometimes) of objective eventful realities of the past.

The current study brings some critical corrections related to the absolute chronology of the KGK-VI communities (humans, not pottery style!), especially regarding the end of this civilization, which cannot be placed around 4200 BC as recently postulated by various scholars, but much later. Why some communities in certain settlements are more resilient than others and the SEB collapse causes remain an open question, which we intend to answer in the future.

Furthermore, the research approach outlined here, must involve in future studies on this topic the identification of potential biases, integration of additional dates, and the correlation with other proxies (e.g., mortality juvenility index, data relative to food production, site and house counts, particular pottery style shift, innovations spreading, climate changes, and paleoenvironment alteration) or other constraints (e.g., possible reservoir effect on ¹⁴C dating on humans and other omnivorous species), which are critical for further developments of SPD analysis of ancient demography in the Balkans and Southeastern Europe in order to increase the accuracy of the data included in the analysis. Additionally, future studies should be extended to incorporate all Chalcolithic signals in the region, thus allowing the correlation of their results to those provided by other studies elsewhere in Europe and the Near East. In this way, it will be possible to better understand the bigger picture of the past demography, human behaviors determinants (necessities, choices, availability, and constraints), and specific social and economic mechanisms that may connect human groups at larger spatial scales in the targeted area.