ADDRESSING THE INTENSITY OF CHANGES IN THE PREHISTORIC POPULATION DYNAMICS: POPULATION GROWTH RATE ESTIMATIONS IN THE CENTRAL BALKANS EARLY NEOLITHIC

Abstract

The intensity of changes in the population dynamics of the Early Neolithic (ca. 6250–5300 cal BC) communities in the Central Balkans was addressed by estimating the growth rate values. The Bayesian approach (Crema and Shoda 2021) of estimating intrinsic growth rates by fitting different models of population growth was applied to radiocarbon dates from the Early Neolithic sites in Serbia. We explored two possible episodes of population growth based on the results of the population dynamics reconstruction using the summed calibrated radiocarbon probability distributions (SPD) method. The results have shown that, within the first episode of growth, the intrinsic growth rate mean values are higher than the estimated continental average (between 1% and 2%). The results indicate a sudden and fast rise in population size, possibly due to the influx of the new population settling in the region at the beginning of the Neolithic. Lower values for the second episode could indicate more gradual population growth due to the mechanisms associated with the Neolithic Demographic Transition and the rise in fertility.

INTRODUCTION

It has been widely attested that the changes that came with establishing the Neolithic way of life happened in all segments of human life–cultural, economic, technological, and biological. A comprehensive theoretical framework for studying all the underlying changes and mechanisms within the Neolithization processes has been formulated by Bocquet-Appel in the form of the Neolithic Demographic Transition theory (NDT in the further text) (Bocquet-Appel 2002; Bocquet-Appel 2008; Bocquet-Appel et al. 2008; Bocquet-Appel and Bar-Yosef 2008; Bocquet-Appel 2011a; Bocquet-Appel 2011b; Bocquet-Appel 2014). The NDT predicts population growth at the beginning of the Neolithic due to the rise in fertility, with the intrinsic growth rate estimated to be between 1% and 2% on a continental level (Bocquet-Appel 2002; Bocquet-Appel and Naji 2006). At one point, mortality rates caught up with fertility, which led to the stabilization of the population size, i.e., the population growth stopped, and the population reflected the state of stationarity, with zero growth rate and constant size through time (Ryder 1975; Bocquet-Appel 2008; Bocquet-Appel 2014).

An important aspect of research and understanding of demographic changes at the beginning of the Neolithic is the reconstruction of population dynamics. Broadly speaking, this includes changes in population size through time as well as the estimation of population growth rates. In the past decade, the method of summed calibrated radiocarbon probability distributions (SPD in the further text) has been widely used as a way of addressing demographic aspects of the Neolithization of Europe (Shennan et al. 2013; Timpson et al. 2014; Porčić et al. 2016; Blagojević et al. 2017; Silva and Vander Linden 2017; Porčić, Blagojević, et al. 2021), but also of other regions worldwide (Lesure 2008; Lesure et al. 2014; Crema et al. 2016). The use of this method enabled the identification of different episodes of growth and decline in population size, enabling researchers to mark distinctive periods of demographic changes while strengthening the main assumptions of the NDT theory. A better understanding of the intensity of these changes could be gained by assessing intrinsic growth rates in the targeted periods.

The aim of this study is to get an insight into the intensity of the already identified demographic changes during the Early Neolithic (ca. 6250–5300 cal BC) in the Central Balkans by estimating the intrinsic growth rates.

Population Dynamics in the Central Balkans Early Neolithic—Previous Studies

Changes in the population dynamics of the Early Neolithic populations in the Central Balkans have been intensively studied in the past several years (Porčić et al. 2016; Blagojević et al. 2017; Porčić et al. 2020; Porčić, Blagojević, et al. 2021; Porčić, Nikolić, et al. 2021). The term “Central Balkans”, in the context of these studies (including this one), geographically covers the territory of modern-day Serbia, while the term “Early Neolithic” refers to the Starčevo culture sequence, covering the time span between ∼6250 cal BC and ∼5300 cal BC (Garašanin and Radovanović 2001; Whittle et al. 2002; Boric 2009; Borić 2011; Porčić et al. 2020). Regardless of the sample used, changes observed in the Central Balkans Early Neolithic population dynamics persisted–they partly resembled the so-called boom-and-bust pattern already observed in different regions in Europe (Shennan et al. 2013; Timpson et al. 2014; Crema et al. 2016; Silva and Vander Linden 2017). This pattern consists of an episode of abrupt growth on the proxy SPD curve at the beginning of the Neolithic, followed by a sudden decline ∼250 years later, after which another growth was identified. This double hump curve has been thoroughly discussed in previous studies, the most recent of which suggested that the first peak could resemble the effect of population migrations from the Near East and the Aegean, while the second peak could represent the effect of in situ fertility growth caused by NDT mechanisms (Porčić, Blagojević, et al. 2021; Blagojević 2022).

The estimation of the Early Neolithic growth rates in the Central Balkans has been done so far by using the Approximate Bayesian Computation (ABC) method on the so-called Probabilistic sample, which consists of 168 dates obtained on animal bones from 27 sites (Porčić et al. 2021). The method was applied for both episodes of growth, by fitting exponential and logistic models of population growth. The resulting growth rates are approximately within the range predicted by the NDT theory; for the first episode of growth, the estimated values were 1.14% for the exponential model and 2.36% for the logistic model (Porčić et al. 2021). Within the second episode of growth, these values were 1.76% when fitting the exponential growth model and 1.92% when fitting the logistic model (Porčić et al. 2021). One problem that emerged, and which we will try to address in this paper, concerns the way time intervals for the analysis were set. More specifically, the chronological span of the model from the previous research had the lower limit set at 6375 cal BC for the first episode of growth (Porčić et al. 2021; Supplementary file 1). In other words, these limits were set too broadly. Since no traces of Early Neolithic settlements in the Central Balkans older than 6250 cal BC have been found so far, the intervals were adjusted according to the available archaeological evidence and results from the palaeodemographic reconstructions.

MATERIALS AND METHODS

In this study, we will use the novel method developed by Crema and Shoda (2021) which is based on growth models with discrete time units that can be fitted through Markov Chain Monte-Carlo (MCMC). The method was originally tested on data regarding the demographic changes related to the period of establishing irrigated rice farming on the island of Kyushu, Japan, during the first millennium BC (Crema and Shoda 2021). This experimental analysis has demonstrated that the method successfully produces real parameters from simulated data in different conditions, even in cases of small sample sizes or different effects of the calibration curve.

The initial sample consists of 331 radiocarbon dates from 53 Early Neolithic sites (Figure 1), some of which are published for the first time in this paper (Supplementary material file 1, Table 2; these specimens were AMS dated at the University of Bristol, School of Chemistry, Bristol Radiocarbon Accelerated Mass Spectrometry–BRAMS), i.e., all currently available Early Neolithic dates from the Starčevo culture sites in Serbia (excluding the Danube Gorges sites due to the specific Neolithization processes characteristic for this region). This sample represents the updated version of the Grand Early Neolithic dates sample (GEN, in the following text) used in previous studies where SPD analyses were applied in population dynamics reconstructions (Porčić et al. 2016; Blagojević et al. 2017; Porčić, Blagojević, et al. 2021; Blagojević 2022). It consists of 235 dated animal remains (70.9% out of the total sample), 54 dated human remains (16.3%), 15 dated charcoal remains (4.5%), and 22 dated plant remains (grains, seeds, fruit stones, and shells; 6.6%). For five samples, this data was unknown.

Figure 1 Map of the sites used in this study: 1. Anište-Bresnica; 2. Autoput E-70, km 521, site 1; 3. Autoput E-70, P2 sever (3); 4. Bakovača-Ostra; 5. Banja-Aranđelovac; 6. Baštine-Obrež; 7. Bataševo; 8. Bezdan-Bački Monoštor; 9. Biserna obala-Nosa; 10. Blagotin; 11. Crnokalačka bara; 12. Crnoklište; 13. Divostin; 14. Donja Branjevina; 15. Drenovac; 16. Golokut-Vizić; 17. Gospođinci-Futog-Klisa I; 18. Gospođinci-Nove Zemlje; 19. Grabovac-Đurića vinogradi; 20. Grivac; 21. Iđoš; 22. Jaričište 1; 23. Kremenilo-Višesava; 24. Kudoš-Šašinci; 25. Lazarev grad-Crkvena građevina; 26. Ludoš-Budžak; 27. Magareći mlin; 28. Međureč-Dunjički šljivari; 29. Miokovci-Crkvine; 30. Motel Slatina; 31. Novi Sad-Gornja šuma; 32. Ornice-Makrešane; 33. Pavlovac-Gumnište; 34. Perlez-Batka C; 35. Pseće brdo-Bečej; 36. Ribnjak-Bečej; 37. Rudna Glava; 38. Rudnik Kosovski; 39. Sajan-Domboš; 40. Sajlovo, site 5; 41. Šalitrena pećina; 42. Selište-Sinjac; 43. Šljunkara na Dumači; 44. Sremski Karlovci-Sonje Marinković; 45. Starčevo-Grad; 46. Staro selo-Idvor; 47. Svinjarička čuka; 48. Topole-Bač; 49. Vinča-Belo brdo; 50. Vinogradi-Bečej; 51. Vršac-At; 52. Zmajevac; 53. Zmajevo-Livnice. The map was produced by T. Blagojević in QGIS 3.32 (www.qgis.org).

The database of Early Neolithic radiocarbon dates in Serbia has been updated with new radiocarbon dates; therefore, a new SPD analysis was performed, mostly to refine the limits of the intervals included in the growth rate estimation analyses. The most recent calibration curve, IntCal20 (Reimer et al. 2020), was incorporated into the analysis. As the resulting pattern on the population curve was not significantly different from the previous results (Porčić et al. 2021), it will not be thoroughly discussed in this study. In this context, the results should be considered a tool for defining intervals and selecting dates to be included in the growth rate estimation analyses (Figure 2).

Figure 2 The results of the SPD analysis on the Grand Early Neolithic sample (N dates: 331, N bins: 108; p<0.001); the dark line represents the empirical curve, the gray area represents 95% confidence intervals based on simulating the SPD curves from the null model (a uniform model that assumes a stationary population, i.e., a constant population size through time), the red areas represent statistically significant growths, and the blue areas represent statistically significant declines on the empirical curve. Intervals used for the growth rate estimation analyses are marked with dark squares for the logistic and with red squares for the exponential growth models for both episodes of growth. The figure was produced by T. Blagojević in the R programming language using the rcarbon package (Bevan and Crema 2018; Crema and Bevan 2020).

Since two distinctive episodes of population growth can be observed on the SPD curve, two separate analyses were performed. In all the analyses, the exponential and logistic growth models were fitted. Based on the shape of the SPD curve as well as single SPD values that fall within the observed statistically significant episodes of growth, an interval for the growth rate estimation was determined. The span of this interval was bounded by the calibrated dates that are included in the first and second episodes of growth, respectively. As a final step, all the dates were included based on their median calibrated values to avoid narrowing the interval too much while staying within the 95% CI limits. Therefore, only subsamples of the total number of 331 dates will be used for the growth rate estimation (Figure 3), although it should be kept in mind that the signals of the episodes of growth are based on the entire sample of radiocarbon dates.

Figure 3 The distribution of calibrated radiocarbon dates (95% confidence intervals) and sites included in the study as a subsample of the Grand Early Neolithic sample for both episodes of growth. The figure was produced by M. Porčić and T. Blagojević in Excel.

Two models of population growth were used in the analysis: the exponential and the logistic growth models. The exponential model was fitted only to the radiocarbon dates falling within the first half of the time interval of an episode of growth detected by the SPD analysis, as it is assumed that the first phase of the growth will be exponential. The formula used for the exponential growth model is defined within the nimbleCarbon package (Crema 2021):

$${\pi _{a - t}} = {{{{(1 + r)}^t}} \over {\sum\nolimits_{t = 0}^{a - b} {{{(1 + r)}^t}} }}$$

where t is time, and a and b are calendar years in BP that define the start and end of the time window of the analysis. The logistic model was fitted to all the dates falling within the interval of an episode of growth detected by the SPD analysis, as the logistic model reflects both the phase of accelerated and decelerating growth (Chamberlain 2006: 21; Porčić 2016: 80). As with the exponential growth model, the formula for the logistic model is also defined within the nimbleCarbon package (Crema 2021):

$${\pi _{a - t}} = {{{1 \over {1 + {{1 - k} \over k}{e^{ - rt}}}}} \over {\sum\nolimits_{t = 0}^{a - b} {{1 \over {1 + {{1 - k} \over k}{e^{ - rt}}}}} }}$$

The first episode of growth covered the time span between ∼6250 cal BC and ∼6000 cal BC, in the case of fitting the logistic growth model, and the time span between ∼6250 cal BC and ∼6125 cal BC when fitting the exponential growth model (Figures 2 and 4). The second episode of growth included dates that fall within the range of ∼5800 cal BC and ∼5650 cal BC, for the logistic model, and from ∼5800 cal BC to ∼5700 cal BC, for the exponential growth model (Figures 2 and 5).

Figure 4 The distribution of calibrated radiocarbon dates (68.2% and 95.4% confidence intervals) from the first episode of growth for logistic and exponential growth models (6250–6000 and 6250–6125 cal BC), plotted on the IntCal20 calibration curve. The figure was produced by T. Blagojević in OxCal 4.4 Online (OxCal).

Figure 5 The distribution of calibrated radiocarbon dates (68.2% and 95.4% confidence intervals) from the second episode of growth for logistic and exponential growth models (5800–5650 and 5800–5700 cal BC), plotted on the IntCal20 calibration curve. The figure was produced by T. Blagojević in OxCal 4.4 Online (OxCal).

The method was implemented in the NimbleCarbon v0.1.0. package, designed specifically for the purpose of estimating growth rates by using different population growth models (Crema 2021).

To account for the research bias, mostly reflected in a large number of dates from some of the individual sites, the binning of dates was performed first. The procedure of binning, which is integrated into the analysis, involves the grouping of the dates within specified time intervals (Shennan et al. 2013; Timpson et al. 2014). Dates are assigned to a bin based on their temporal distance, which is set by a researcher. In this case, it was set at 150 years. In other words, if the distance between two dates is up to 150 years, they would belong to the same bin. One date was randomly chosen from each bin (provided that it is within the temporal limits of the analysis). This step affected the effective sample size within the defined intervals by lowering the number of dates included. All the dates that were used in the analysis are given in Table 2, in Supplementary material file 1, while the number of dates that entered the final sample is given in Table 1. The code for the analysis in R is given in the Supplementary material file 1. This procedure was repeated multiple times (10 times) for both episodes and for both models to check for sampling effects. The results of all iterations are given in Table 1 in Supplementary material file 1. Supplementary material file 2 is an Excel file that contains almost all the dates that were used in the analysis, arranged in a way to be directly loaded in R using the code provided in Supplementary material file 1, and for the analysis to be reproduced. It should be stressed that nine dates from the original analysis were not included in this file since their Uncal BP values have not been published yet. However, given the sample size, it can be assumed, with a high degree of certainty, that their omission won’t affect the overall pattern.

Table 1 The results of the growth rate estimations for different episodes of growth.

The priors for parameters of the models were also defined: r (growth rate)—for both exponential and logistic models; and k (initial population size, which, in this case, represents the proportion of carrying capacity, K), only for the logistic growth model. The parameter values were defined within the Bayesian framework, which implies confronting the empirical data with prior probability distributions of parameter values and generating posterior distributions. The prior for parameter r (growth rate) was modeled as a uniform distribution at the interval between 0.0001 and 0.06, while the prior for parameter k (the ratio of initial population size and carrying capacity in the logistic model) was modeled as a uniform distribution at the interval between 0.0002824891 and 0.01129956. If we assume that the carrying capacity in the Early Neolithic Balkans was one person per square kilometer (Dennell and Webley 1975; Müller 2006; Lemmen 2013), given the area of the Republic of Serbia, which is 88499km², the lower limit value for k corresponds to the assumption that there were 25 people initially, whereas the upper limit value is based on the assumption that the initial population size was 1000 people. These estimates are little more than guesses, but the prior distribution for k should be wide enough to encompass all possibilities.

A goodness-of-fit for the models was evaluated by looking at Pearson’s correlation coefficient, which measures the strength of a linear correlation between modeled and empirical SPD values (Crema 2021; Crema and Shoda 2021).

Potential Study Limitations

Since the growth rate estimation analysis was done based on the results of the SPD analysis, several limitations of this method should be addressed in order to ensure the interpretation of the results is as objective and robust as possible.

Biases that can affect the results of the SPD method are sampling bias, calibration curve effects, and taphonomic bias (Williams 2012; Shennan et al. 2013; Contreras and Meadows 2014; Timpson et al. 2014; Heaton 2022).

Sampling bias refers both to the sample size and to decisions related to taking samples from different sites or specific contexts within them, which depend on different research questions. They can be directed towards the reconstruction of the duration of the entire settlement or attempts to date some particular events. Furthermore, specific contexts may be dated multiple times by different researchers, resulting in a significantly greater number of dates compared to other contexts from the same period.

In order to reduce the effect of sampling bias, the SPD method uses a binning procedure. This procedure is particularly important in cases where there are a large number of dates from a certain site, and it implies the chronological grouping of all dates originating from it. Bins are constituted by those dates whose values are chronologically close. More specifically, dates are distributed into new bins when there is a certain chronological difference between them, as defined by the researcher. Most often, it is about 100 or 200 years (Shennan et al. 2013; Timpson et al. 2014). In this research, h=150 years. In this way, “phases” are formed within the site, which consist of dates between which the difference is not greater than 150 years. This ensures that each resulting “phase” has equal weight when generating the final population curve. More specifically, this procedure allows for the reduction of bias within the already existing sample. However, it cannot affect the biases that result from decisions about the selection of sites and contexts for dating. After the binning, the effective sample size is smaller in relation to the total number of dates because now the basic unit of measure is not each date individually but each of the bins.

The calibration curve effect primarily refers to the dependence of individual calibrated values, that is, the shape of their probability distributions, on the shape of the calibration curve at a given moment in time. More precisely, how steep or flat the calibration curve is at a particular time interval will directly affect the probability distribution of individual radiocarbon dates (Williams 2012; Weninger et al. 2015; Crema and Bevan 2020; Heaton 2022). The method used to produce the SPD curve used in this paper takes into account the effects of the calibration curve (Crema and Bevan 2020). This, however, does not mean that the method smooths out the effects of the calibration curve completely. The direct consequence of the plateau observed within the time frame this research focuses on is wider confidence intervals and, hence, less precise results. Therefore, the interpretation of the shape of the curve should be done in the context of other indicators of possible fluctuations, both in population dynamics and in the calibration curve.

Taphonomic bias refers to the influence of taphonomic processes on the archaeological record in its broadest sense, leading to a tendency to lose information over time, typically following an exponential trend. In earlier studies for the area of the Central Balkans (Porčić et al. 2016; Blagojević et al. 2017), the so-called universal taphonomic model (Surovell et al. 2009) was used. However, since the shape of the taphonomic curve within the time interval investigated in this paper indicates that taphonomic processes did not have a significant impact on the obtained results, it was decided not to use a universal taphonomic model as a null model but rather a uniform model that assumes a constant population size through time (i.e., a stationary population) (Crema and Bevan 2020). This certainly does not imply that taphonomic processes do not have a significant role in the formation of the archaeological record, but only that, within the framework of this research, it is possible to minimize their influence.

RESULTS

The results of the growth rate estimations are given in Table 1. The r values in the table represent an estimated growth rate mean value, shown in percentages, with 95% confidence intervals.

The mean values of the estimated growth rates for the first episode of growth are 3.25% for the logistic model and 3.22% when fitting the exponential model (Table 1, Figure 6). The absolute fit of the model was obtained by generating the SPD values from the fitted posterior models, which were then visually compared to the empirical SPD values as a type of posterior predictive check (Crema 2021) (Figure 7). In this case, the curve indicates a good fit of both models. Pearson’s correlation coefficient is high in both cases, indicating a strong correlation (Table 1).

Figure 6 Distribution of growth rate values (r) for (a) logistic and (b) exponential models for the first episode of growth (6250–6000 and 6250–6125 cal BC) for the GEN sample. The figure was produced by T. Blagojević in the R programming language, using the nimbleCarbon package (Crema 2021).

Figure 7 The goodness-of-fit for (a) logistic and (b) exponential growth models for the first episode of growth (6250–6000 and 6250–6125 cal BC) for the GEN sample. The black line represents the empirical SPD values, and the gray shaded area indicates the range of potential SPD curves generated by the models with parameter values sampled from posterior distributions and with equal sample size as the empirical data set. The figure was produced by T. Blagojević in the R programming language, using the nimbleCarbon package (Crema 2021).

For the second episode of growth, the mean values of the growth rate for the logistic model are 1.89% and 2.36% for the exponential model (Table 1, Figure 8). The curves that indicate the absolute fit show a good fit of both models with a short interval on which the curve exceeds the upper limits of the model in the case of fitting the exponential model (marked in red) (Figure 9). This is confirmed by Pearson’s correlation coefficient values, which are lower in the case of the exponential model and indicate a moderate correlation (Table 1).

Figure 8 Distribution of growth rate values (r) for (a) logistic and (b) exponential growth models for the second episode of growth (5800–5650 and 5800–5700 cal BC) for the GEN sample. The figure was produced by T. Blagojević in the R programming language, using the nimbleCarbon package (Crema 2021).

Figure 9 The goodness-of-fit for (a) logistic and (b) exponential growth models for the second episode of growth (5800–5650 and 5800–5700 cal BC) for the GEN sample. The black line represents the empirical SPD values, and the gray shaded area indicates the range of potential SPD curves generated by the models with parameter values sampled from posterior distributions and with an equal sample size as the empirical data set. The figure was produced by T. Blagojević in the R programming language, using the nimbleCarbon package (Crema 2021).

DISCUSSION AND CONCLUSION

Several different studies that used the SPD method for population dynamics reconstruction in the Central Balkans Early Neolithic produced a general pattern that persisted regardless of the samples used (Porčić et al. 2016; Porčić, Blagojević, et al. 2021; Blagojević 2022). This pattern includes an abrupt and significant population growth at the beginning of the Early Neolithic that lasted for ∼250 years, followed by a sudden drop in the population curve. The drop lasted for about 200 years, after which a new significant growth was detected (Porčić et al. 2016; Porčić et al. 2021; Blagojević 2022). Intrinsic growth rate values were estimated for two episodes of growth.

For the first episode of growth, the intrinsic growth rate mean values are considerably higher than the estimated average for the NDT (between 1% and 2%) (Bocquet-Appel 2002), both when fitting the logistic and the exponential models (Table 1). As for the second episode of growth, the estimated values broadly correspond to the NDT expectations, even though they are slightly higher than the estimated upper limit (Table 1).

In the previous study (Porčić et al. 2021), the estimated growth rate values were lower for both episodes of growth. A significant difference between the previous and the results obtained in this study could be primarily explained by differences in the method. These include both the process of defining the interval and the methods used.

The results of this study indicate a very rapid and abrupt increase in population size, with high growth rates, at the beginning of the Neolithic, which can only partially be explained by an increase in fertility alone. One of the possible explanations for the higher end of the spectrum of the estimated growth rates, would be that these results represent the signal of a mixed effect of the influx of a new population to the previously unsettled area during a short time interval, as well as local population growth.

Lower growth rate values estimated for the second episode of growth might be considered in light of previous interpretations that this episode of growth actually represents the signal of in situ population growth due to the rise in fertility caused by the mechanisms associated with the Neolithic demographic transition (Dennell and Webley 1975; Lemmen 2013; Porčić et al. 2021).

Another important implication of this study relates to how the sample size affects the results. Based on this study alone, it cannot be argued whether there is a minimal sample size required for the analysis to be robust. However, it could be seen that the least fit of the models, as seen through Pearson’s correlation coefficient values, was obtained with the smallest sample (the second episode of growth, an exponential model) and that this difference is high once the number of dates goes below 10.

As with any other estimations regarding the demographic changes and their dynamics during the Early Neolithic, it should be again emphasized that they represent a valuable tool in detecting important patterns. Nevertheless, these estimations should primarily be considered as an insight into particular changes. This could represent an important frame for constructing future research questions that should be assessed only through a multi-scale approach that would include multiple (available) archaeological indicators.