The Neolithic Demographic Transition in the Central Balkans: population dynamics reconstruction based on new radiocarbon evidence

In this paper, we test the hypothesis of the Neolithic Demographic Transition in the Central Balkan Early Neolithic (6250–5300 BC) by applying the method of summed calibrated probability distributions to the set of more than 200 new radiocarbon dates from Serbia. The results suggest that there was an increase in population size after the first farmers arrived to the study area around 6250 BC. This increase lasted for approximately 250 years and was followed by a decrease in the population size proxy after 6000 BC, reaching its minimum around 5800 BC. This was followed by another episode of growth until 5600 BC when population size proxy rapidly declined, reaching the minimum again around 5500 BC. The reconstructed intrinsic growth rate value indicates that the first episode of growth might have been fuelled both by high fertility and migrations, potentially related to the effects of the 8.2 ky event. The second episode of population growth after 5800 BC was probably owing to the high fertility alone. It remains unclear what caused the episodes of population decrease. This article is part of the theme issue ‘Cross-disciplinary approaches to prehistoric demography'.


Introduction
Several major theories and models have been proposed to explain the Neolithic expansion across Europe. For this reason, palaeodemographic research on a regional scale is needed to explore the process at higher spatial and temporal resolutions, as we now have the methods and growing body of data to test the models and theories. In the past decade, probably the most popular method for the reconstruction of population dynamics was summed calibrated probability distributions (SCPD) of radiocarbon dates. This method had been originally proposed by Rick [1]. The two basic assumptions for this method are that (i) there is a correlation between population size and quantity of material remains and (ii) the sample of material remains used for dating is representative (it approximates a random sample). The implication is that the frequency distribution of absolute dates should reflect changes in relative population size ( population dynamics) over a time interval. As radiocarbon dates are not point estimates (e.g. a year) but probability distributions, the dates cannot be simply counted, but their calibrated probability distributions are summed in order to derive a distribution of probabilities associated with each calendar year in the interval. This summed probability curve is interpreted as a relative population size proxy. The method has been criticized on many grounds-the assumptions were questioned, especially the representativeness of the radiocarbon dates owing to research bias, and it was pointed out that many non-demographic factors such as changes in settlement patterns, taphonomic loss, sampling error and the shape of the calibration curve might influence the shape of the SCPD curve [2][3][4]. The problems of between-site research bias in the number of radiocarbon dates, taphonomic loss, sampling error and influence of the calibration curve were mostly solved or significantly reduced in the modern formulation of the SCPD method by Shennan et al. [5] and Timpson et al. [6]. As for the influence of other factors, there is no denying their reality, but rather than rejecting the entire approach, the way forward is to evaluate their influence empirically on a case-by-case basis [7].
In parts of Europe where a large body of dates come from a commercially driven archaeology and where the quantity of radiocarbon dates is large in general (e.g. northern, western and parts of central Europe), the assumptions of the SCPD method are met or approximated to a reasonable degree (but see Torfing [3][4] for a different opinion). One of the earliest applications of the modern version of the SCPD method was to reconstruct the population dynamics of the Neolithic populations in northwestern and central Europe-the results revealed the boom-and-bust pattern across different regions [5,6]. The rapid increase of population size at the beginning of the Neolithic was in line with the predictions of the major theories of the Neolithic expansion, but the subsequent population crash was not. Shennan et al. [5] speculated that the 'bust' episode could be interpreted in Malthusian terms as a consequence of soil depletion and erosion owing to unsustainable farming practices and overpopulation. More recently, it has been shown that the cause of the population crash in northwestern Europe was owing to the agricultural crisis [8]. These studies have shown the importance and potential of regional investigations for understanding the process of the Neolithic expansion. In this paper, we move the focus to southeastern Europe, specifically to the Central Balkans.
Until recently [9][10][11], there was practically no information about the population dynamics of the Neolithic populations closer to the origins of the Neolithic expansion in Europethe Central Balkans. The Early Neolithic communities from the Aegean reached the Central Balkans by approximately 6250 BC according to the new radiocarbon dates published in this study (see also [12]). Traditionally, this Early Neolithic has been labelled as the Starcěvo culture, which is a part of the wider Early Neolithic cultural Starcěvo-Körös-Criş complex, lasting from approximately 6250 to approximately 5300 BC [13][14][15][16]. The Central Balkan region was an important station for the Neolithic expansion into Europe, as farmers coming from the southeast in the second half of the 7th millennium BC faced somewhat different ecological conditions, which might have required adaptations that enabled the further spread of the Neolithic [17].
Previous reconstructions of the Central Balkan population dynamics were based on the application of the SCPD method to a small sample of published radiocarbon dates [9]. The results revealed two episodes of rapid population growth (between 6200 and 6000 BC and between 5800 and 5600 BC), separated by an episode of substantial population decrease between 6000 and 5800 BC.
In this paper, we present the results of palaeodemographic research aimed at reconstructing the population dynamics of the first farmers of the Central Balkans, using the SCPD method applied to the sample of over 200 new radiocarbon dates. Unlike previous applications of the method with datasets collated from the literature, our sample was designed and collected specifically for the application of the SCPD method. The fact that we were able to take a random and representative sample effectively eliminates any research bias towards or against particular sites or periods, and thus removes one of the major objections to the SCPD method (for critical reviews of the SCPD method see [2,3]). Our proximate aim is to reconstruct population dynamics of the first farmers in the Central Balkans between 6250 and 5300 BC, and our ultimate aim is to discuss these results in the light of the major theories and hypotheses about the Neolithic expansion in Europe, primarily the Neolithic Demographic Transition theory.

Major theories and models of the Demic Neolithic expansion
The Wave of Advance (WoA) model formulated by Ammerman and Cavalli-Sforza [18] is a theoretical explanation for the spread of the Neolithic from the Near East into and through Europe. It assumes that local Neolithic populations were growing according to the logistic growth curve. The logistic model of population growth is a model in which the population starts growing with increasing speed, but once it reaches a certain point, the growth starts to decelerate and it finally stops once the population reaches the carrying capacity [18,19]. In parallel with this process, the model predicts that people migrate to the neighbouring regions, provided that they are not populated up to their carrying capacities. The model assumes no particular direction for the migration, as all directions in the immediate neighbourhood of the settlement are equally probable a priori, but the aggregate effect is a directional movement of the farming front away from the origin and towards demographically empty space as the immediate surroundings of each settlement get filled up by farmers. If this model were true, we would expect to see the SCPD curve from the Central Balkan Early Neolithic following the logistic model. The theory of the Neolithic Demographic Transition (NDT) as formulated by Bocquet-Appel [20] explains why Neolithic populations would experience population growth in the first place. It has been known for a long time that the Neolithic way of life was linked with rapid population growth, but researchers disagreed on whether this growth was caused by the increased fertility (e.g. [21]) or decreased mortality (e.g. [22,23]). The theory of the NDT suggests that the increase in population size was caused by the increase in the total fertility rate [20,24,25]. The core of the theory is the Relative Metabolic Load Model, in which the number of children depends on the amount of available energy determined by the general conditions of life [20]. According to this theory, when people adopted a farming way of life in the Neolithic, the reduction of mobility and the reliance on a diet rich in carbohydrates cut down energy expenses and increased energy inputs for females, resulting in a higher amount of energy available for reproduction. As fertility increased, farming populations started to grow at relatively high rates. This is the first phase of the NDT and the signal of increased fertility has been empirically detected in numerous studies [26][27][28]. The second phase of the NDT started when the high population density in settlements, close contact with animals, and a diet rich in calories but poor in nutrients resulted in health deterioration owing to infectious diseases and nutritional stress, ultimately leading to an increase of mortality [20]. At one point, the mortality rate became equal to the birth rate and the further population growth was stopped. The theory of the NDT provides a clear prediction for the Neolithic population dynamics. It should also follow the logistic model, as this model captures the rapid increase owing to high fertility, and a gradual deceleration of growth when mortality equalized with fertility.
There are also other hypotheses about the spread of the Neolithic that are consistent with genetic evidence, but are different from the WoA and NDT models. Recently, Silva & Vander Linden [11] proposed the travelling wave-front (TWF) model as an explanation of the boom-and-bust patterns in the European Neolithic. This model assumes that once the Neolithic population had arrived in some region by largescale migration, this would be reflected as a sharp increase in the SCPD proxy, whereas a sudden decrease of the SCPD curve is interpreted as, 'due to outgoing migrants, resuming their spread into a new region' [11, p.6]. Although it also includes a demographic wave, the TWF model is different from the WoA model. In the WoA model, the population in each region does not decrease after a population wavefront moves to a new region, whereas in the TWF model, the population of the region decreases after the wavefront moves onwards. Additionally, the wave signal in the TWF model is mostly a result of large-scale migrations rather than in situ growth owing to increased fertility.
Finally, we mention the hypothesis proposed by Weninger et al. [29], which is not a model or theory, but a historical scenario for the spread of the Neolithic. The main idea is that changes in climate associated with the 8.2 ky event caused the migration of the farming communities from Anatolia into Europe [29]. The 8.2 ky represents an episode of an abrupt climate change in the Holocene, resulting in decreased temperatures between approximately 8200 and approximately 8000 ky calBP (approx. 6250-6050 BC) with potential disrupting effects on the terrestrial ecosystems and the Neolithic farming systems [29,30].

Materials and methods
The study area is the Central Balkans (figure 1) and the temporal range of the study is the period from approximately 6250 BC to approximately 5000 BC. There are three partially overlapping sets of radiocarbon dates from the Republic of Serbia (core area of the Central Balkans) used in this study (for details of sampling see electronic supplementary material, file S1; all radiocarbon dates used in this study are in electronic supplementary material, file S2): (1) The probabilistic sample consists of 168 accelerator mass spectrometry (AMS) radiocarbon dates from 27 sites collected by the European Research Council BIRTH project (see Acknowledgements) with an idea to approximate random sampling. This sample is collected exclusively for the purpose of SCPD palaeodemographic reconstruction in order to eliminate research bias towards particular sites or periods within the Early Neolithic sequence. (2) The grand Early Neolithic sample consists of 1. the probabilistic sample, 2. a sample of radiocarbon dates produced by the BIRTH project for purposes other than SCPD reconstruction, and 3. previously published Early Neolithic dates collated in [9]. This sample consists of 299 dates from 52 Early Neolithic sites in Serbia. (3) The human remains sample. As all dates from the probabilistic sample are based on animal bones, we decided to use the sample of radiocarbon dates made only on human bone material as partially independent ( partially because specimens in both samples come from the overlapping set of sites) of the probabilistic sample. There are 45 radiocarbon dates on human remains from 20 sites in this sample

(a) Analysis
In order to reconstruct population dynamics, the SCPD method was applied to each of the three radiocarbon datasets. We used the standard SCPD method as formulated in [5,6], supplemented by the extension that tests for the statistical significance of pointto-point differences on the SCPD population proxy curve, as published in [31]. The empirical SCPD curve is compared to the 95% confidence intervals based on the set of SCPD curves generated by sampling from the null model. The null model assumes that the underlying population was stationary (constant population size) throughout the time interval of interest. When the empirical SCPD is above or below the 95% confidence intervals, there is a statistically significant growth or decline of the population relative to the null model. The methods were implemented in the rcarbon package [32] for R [33] (for technical details of the analysis see electronic supplementary material, files S1 and S4; data are stored in the electronic supplementary material, file S2). In order to estimate intrinsic growth rates based on the SCPD curve, we used the Approximate Bayesian Computation (ABC) method [34,35]. The standard regression techniques are not appropriate in this case as the sample size of the SCPD curve (one point per year) would be artificially inflated. In addition to this problem, the calibration curve may strongly distort the SCPD curve in time intervals where the curve is flat or wiggly. For this reason, the regression would not be able to give accurate parameter estimates.
We assumed that the population was growing according to an either logistic or exponential model. The aim of the ABC analysis is to estimate the growth rate for each model. This is accomplished by simulating a large number of population dynamics scenarios that differ in the intrinsic growth rate values (sampled from a prior distribution of possible values). For each simulated scenario, a SCPD curve is produced and compared to the empirical SCPC curve for the goodness of fit. Intrinsic growth rate values that generate simulated SCPD curves with a good match to the empirical curve are used to generate a posterior distribution of the growth rate estimates. For technical details of the analysis, see electronic supplementary material, files S1 and S4. This analysis is performed only on the SCPD from the probabilistic sample based on calibration with normalization.

(a) Population dynamics reconstruction based on the summed calibrated probability distributions
The results of the SCPD analysis for the probabilistic sample are shown in figure 2. Regardless of the presence or absence of normalization [36], the resulting pattern is the same. The curve increases significantly from approximately 6300 BC, reaches its peak around 6050 BC and then significantly decreases, reaching its minimum around 5800 BC. Afterwards, the curve starts to increase again, reaching a peak around 5600BC and again plummeting beneath the lower 95% confidence interval limit between approximately 5500-5400 BC. This is followed by another increase of the curve after 5400 BC reflecting the population increase associated with the start of the Late Neolithic. Therefore, we have two major episodes of statistically significant growth, from approximately 6250 to approximately 6000 BC and from approximately 5800 to approximately 5600 BC, respectively, royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 376: 20190712 and two major episodes of decline, from approximately 6000 to approximately 5800 BC, and from approximately 5600 to approximately 5400 BC, respectively. The global significance is below the 0.001 level, indicating that these deviations from the null model are highly unlikely to be owing to chance. The almost identical pattern is found for the grand Early Neolithic sample and human remains sample, with global significance in each case being lower than 0.005. For these two samples, the trough in the curve around 5800 BC does not go under the lower limit of the 95% confidence intervals, but the specific statistical test [31] for changes in the shape of the SCPD curve suggests that differences between the highest point around 6000 and lowest point around 5800 BC are statistically significant at the 0.005 level.

(b) Estimating growth rate values
As there are two major episodes of growth of the SCPD curve, the estimation was performed for each episode separately (see electronic supplementary material, file S1). For the first episode of growth (between approximately 6250 and approximately 6050 BC), the estimated value of the intrinsic growth rate for the exponential model is 1.14% (95% credible interval: 0.61-1.71%) and 2.36% (the 95% credible interval: 0.9-3.97%) for the logistic model. For the second episode of growth (between approximately 5800 and approximately 5600 BC), the estimated value of the intrinsic growth rate for the exponential model is 1.76% (95% credible interval: 0.98-2.57%) and 1.92% (the 95% credible interval: 0.89-3.1%) for the logistic model.

Discussion
The results suggest that there was an increase in population immediately following the introduction of the Neolithic into the Central Balkans, as predicted by all models. This increase lasted for around 250 years and was followed by another episode of the Early Neolithic population growth, culminating around 5600 BC.
Estimates of the NDT growth rates by Bocquet-Appel [27] are usually between 1 and 2%. Estimates in this study for the exponential model are within the range of the NDT values. For the logistic model, the estimated values are closer to the higher end of the spectrum of the recorded human growth rates [18]. However, in cases where population moves through space and time, the resulting aggregate curve that is based on mixing data from successively colonized regions may have different properties from local growth curves (see electronic supplementary material, file S1). When population proxy data from different regions are combined, the resulting space-averaged curve is less steep than individual local growth curves. The implication is that the local population growth rates are actually much higher than growth rate estimates based on the space-averaged aggregate curve. This would mean that the local growth rates for the first episode of growth when the farming frontier was moving to the north were even higher than the reported estimates, on the very upper (or over) the limit for the values expected for the NDT. A recurring pattern of the Neolithic population booms followed by busts was found in many regions of Western and Central Europe [5,6]. The TWF interpretation is logically possible for the episode of population decrease between approximately 6000 BC and approximately 5800 BC-in this interval the SCPD curve in Central Balkans goes down, while in the Great Hungarian Plain, it goes up [10], which is consistent with the TWF model. Therefore, the Central Balkan evidence is consistent with this model, at least for the first episode of population boom-and-bust. If this model is accepted for the explanation of the observed pattern, the first episode of growth would reflect the arrival of the Neolithic travelling wave-the observed increase of population between 6250 and 6000 BC would be owing to the large-scale migration to the region. Most of this population would migrate after approximately 6000 BC, leaving the Central Balkans at the low population level. The second episode of growth, starting after approximately 5800 BC, would represent the in situ population growth owing to the NDT.
The boom-and-bust pattern from 6250 to 5800 BC is consistent with the TWF model, but the results are also consistent with the WoA model-the first 200-250 years of the Neolithic expansion occurred consistently with the WoA, but the unknown exogenous factor caused the population decrease after 6000 BC. Currently, there is no evidence for epidemics or intensive violence around 6000 BC that might indicate a collapse caused by overpopulation, although this may more reflect the absence of evidence in general, rather than evidence of absence. An additional problem with the TWF model is that it is mostly descriptive, as it is based on empirical observations. There is nothing in this model that tells us what is the endogenous engine of the migrations.
The 8.2 ky event hypothesis would seem to be a good candidate for the cause of the chain migrations. If the historical Little Ice Age is taken as a relevant analogy for the 8.2 ky event, then we should expect many adverse effects on the farming system [37]. The population growth in the Central Balkans was counterintuitively most intensive between approximately 6250 and approximately 6050 BC (or approximately 6000 BC, depending on the sample), coinciding for the most part with the 8.2 ky event. If we follow the 8.2 ky migration hypothesis [29], the event caused a massive migration of farming communities from the Aegean towards the north. This influx of migrants pushed the population size of the colozined regions over the carrying capacity in a very short time, as the natural growth of the Neolithic population owing to increased fertility would be augmented by the migration component. The resultant growth rate would be very high, which is consistent with the results of this study. There is evidence that expansion rates increased during the 8.2 ky event [38]. The overshooting of the carrying capacity would trigger further migrations and might have led to population collapse owing to resource depletion or conflict. But little is known about the actual effects of the 8.2 ky event on the cultural systems in Europe. A recent study has shown that it did not systematically affect the Neolithic and Chalcholithic communities in the Near East [39]. In Greece, there are indications that the 8.2 ky event had taphonomic impact in the form of stratigraphic disturbances on Neolithic sites possibly associated with increased intensity of summer thunderstorms, but as the authors of the study stated, more research is needed for firm conclusions [40]. Furthermore, if there was a massive migration from the Aegean triggered by the environmental effects of the 8.2 ky event, this should be reflected in the population proxy as a major decrease of population. If we look at the SCPD curve from Greece, we can indeed see a small trough between 6200 and 6000 BC [41], but without the specific test it is not possible to tell whether this decrease is statistically significant. The link between the 8.2 ky event and demographic aspects of the Neolithic expansion in the Central Balkans is worthy of further investigation.
Finally, a non-demographic explanation is possible for the observed pattern of SCPD decrease after 6000 BC. If residential mobility was higher before approximately 6000 BC and was reduced afterwards, this could produce the pattern of boom-and-bust in the SCPD curve even if the population size was equal. The same effect would be observed if the population agglomerated in a smaller number of larger settlements. Currently, there is no convincing evidence for the change in residential mobility or settlement size after 6000 BC (electronic supplementary material, file S1), but this possibility needs to be investigated further.

Conclusion
In this study, we reconstructed the population dynamics in the Central Balkans during the first thousand years of the Neolithic. The population size constantly and rapidly increased in the first 200-250 years following the arrival of agriculturalists. Estimated growth rates indicate that this increase was on the very edge or over the edge of what can be expected from the increase of fertility owing to the NDT, which indicates that this initial growth might have been fuelled both by high fertility and migration. The potential migrations might be linked to the effects of the abrupt cooling associated with the 8.2 ky event. The first episode of population growth stopped around 6000 BC when the population size proxy began to decrease significantly, reaching its minimum around 5800 BC. It is unknown what caused this decrease. The population quickly recovered by 5600 BC. This second episode of population growth was probably owing to fertility alone. After this second episode of growth, another decrease followed, again with an unknown cause. The research challenge for the future is to explore whether this cyclical pattern of booms and busts is caused by some deeper principle of long-term (pre)history, or the observed boom-and-bust cycles are historical contingencies, each caused by different factors.
Ethics. The permit for the export of the Neolithic bone samples from the Republic of Serbia for the purpose of radiocarbon dating is provided as a part of the documentation (Annex 5) for the BIRTH project (Grant Agreement 640557).
Funding. This study was funded by the European Research Council within the European Union's Horizon 2020 research and innovation programme (Grant Agreement no. 640557).