Modelling the initial expansion of the Neolithic out of Anatolia

Computer-based simulations of the Neolithic expansion in Eurasia using the time-space distribution of 14C dates have consistently highlighted a gradient from the Near East to the British Isles (Gkiasta et al. 2003; Pinhasi et al. 2005; Bocquet-Appel et al. 2009; Fort et al. 2012). The underlying assumption that agriculture swept across Europe following the advance of a pioneer front has (if anything) comforted Childe’s and other diffusionists’ accounts of a migration of early culture based on similarities in the pottery and other material remains (Childe 1925; 1950; Elliot Smith 1915[1929]). Clark is widely credited with the first explicit use of radiocarbon dates for modelling Neolithic expansion (Clark 1965). What has changed since Clark, as other authors have pointed out, is not so much the scope as the resolution of the model, which has improved dramatically thanks to the widespread use of 14C dating (Bocquet-Appel et al. 2009.807).


Introduction
Computer-based simulations of the Neolithic expansion in Eurasia using the time-space distribution of 14 C dates have consistently highlighted a gradient from the Near East to the British Isles (Gkiasta et al. 2003;Pinhasi et al. 2005;Bocquet-Appel et al. 2009;Fort et al. 2012). The underlying assumption that agriculture swept across Europe following the advance of a pioneer front has (if anything) comforted Childe's and other diffusionists' accounts of a migration of early culture based on similarities in the pottery and other material remains (Childe 1925;1950;Elliot Smith 1915[1929). Clark is widely credited with the first explicit use of radiocarbon dates for modelling Neolithic expansion (Clark 1965). What has changed since Clark, as other authors have pointed out, is not so much the scope as the resolution of the model, which has improved dramatically thanks to the widespread use of 14 C dating (Bocquet-Appel et al. 2009.807).
The sheer number of published radiocarbon dates is such that we advocate moving a step further, by drawing a regional simulation of the Neolithic dispersal, this time within a moderately small section of Eurasia (c. 1 000 000km 2 ), spanning from the Central Anatolian Plateau in the East to Thessaly in the West, and the Balkan Range in the North (Fig.  1). Sites in Northern Bulgaria and Serbia fall outside the scope of this paper and will not be considered further here. In the article, we use empirical Bayesian kriging to interpolate the advance of the Neolithic from the Anatolian heartland to Southeast Europe. The model relies upon a comprehensive dataset of 71 sites and 1162 uniformly recalibrated dates, falling within the interval 9000-5500 calBC at 2σ. Unlike other simulations of the Neolithic, which use the oldest observed 14 C date(s) as a proxy for the advance of a pioneer front (e.g., Pinhasi et al. 2005;Bocquet-Appel et al. 2009;Fort et al. 2012), our simulation draws upon modelled dates, statistically constrained by prior information using Bayesian clustering. It goes without saying that this approach is feasible only with a small sample of sites, over which strict quality control can be maintained.
The central question being asked of the data is whether the spread of the Neolithic out of Anatolia was a linear process, or whether it consisted instead of standstills, punctuated by rapid advances. What is at stake is the potential identification of so-called farming 'frontiers' within the study region similar to the ones identified in the Great Hungarian Plain (Whittle 1996;Zvelebil, Lillie 2000), the southern Adriatic coast (Forenbaher, Miracle 2006), the circum-Baltic region (Whittle 1996;Zvelebil 1998; and the Low Countries (Louwe Kooijmans 2007). The traditional view, held by Ammerman and Caval-li-Sforza, is that farming expanded across Europe at a steady pace of approx. 1 km/year (Ammerman, Cavalli-Sforza 1971;1984.61, 135). This estimate, which has been upheld in recent literature (Pinhasi et al. 2005), is at odds with the archaeological picture outlined above and recent demographic work, which suggests an expansion in 'booms and busts' (Shennan, Edinborough 2007;Shennan et al. 2013). The latter pattern of spread is usually captured under the concept of 'arrhythmic' expansion (Guilaine 2000).
By challenging the linear narrative of farming expansion within the study region, we hope to contribute to a growing body of literature which highlights the crucial role of Anatolia not just as a land bridge, but also as an independent centre of neolithisation (Özdogan, Basgelen 1999;Özdogan et al. 2012;Thissen 2000;Gérard, Thissen 2002;Lichter 2005;Gatsov, Schwarzberg 2006;Krauß 2011;Baird 2012;Çilingiroglu 2012). One of the key issues emerging over the years has been the distinction of two Neolithic traditions, one in Central Ana- tolia, running broadly concurrent with Pre-Pottery Neolithic B societies in the Near East, and the other in Western Anatolia, coinciding or shortly pre-dating the widespread adoption of pottery in the Northern Levant (Schoop 2005;Baird 2012;Düring 2013). As this study demonstrates, the advent of farming in Western Anatolia was delayed by up to 2000 calibrated years and this lag in the dating needs to be properly accounted for in future.

Dataset and methods
A geostatistical (kriging) method was used to interpolate the spatiotemporal advance of the Neolithic from a set of known values. The first step was to obtain the known values from the sample data -a georeferenced dataset of 1162 calibrated radiocarbon dates from 71 sites (Electronic Supplementary Material 1). This number excludes duplicate entries and dates that fall outside the range 9000-5500 calBC at 2σ. For the period under review, 1057 dates were ascribed to Neolithic and Early Chalcolithic levels, 99 to Epipalaeolithic and Mesolithic levels; 6 came from mixed layers or could not be ascribed to a period in particular. A Bayesian model was built for each site where it is possible by using median estimators of phase boundaries in OxCal 4.2. Two versions of the kriging, one including virtually all modelled dates, regardless of quality, the other based on a strictly audited sample, were constructed. In turn, the intensity of the Neolithisation process was evaluated through summed probability distributions of calibrated radiocarbon dates.

C data collection, calibration and quality control
The radiocarbon database on which this study relies was collated from published literature and existing datasets, including the CalPal-database (Weninger 2014), the CONTEXT database (Böhner, Schyle 2008) and the CANeW (Reingruber, Thissen 2005;Thissen 2006;Gérard, Thissen 2002). Dates were uniformly recalibrated using the IntCal13 atmospheric curve (Reimer et al. 2013) in Oxcal 4.2 (Bronk Ramsey 2013). The consistency of the database was checked for out-of-scope and duplicate entries. In attributing sites or phases to the 'Neolithic', we followed the assessment of the excavators, cross-checking (where possible) the validity of this attribution, based on such criteria as the adoption of food production, e.g., domestic plants and/or animals (Childe 1936). On this basis, three of the 71 sites surveyed did not return any 'Neolithic' dates and were not processed any further. Subsequently, two different approaches were pursued. The first one involved limited pre-sorting, excluding only those radiocarbon determinations reported as problematic by the laboratories. The advantage of this method is that virtually all 14 C dates, regardless of quality, could be included in the model, thus pre-empting biases regarding the way in which the selection was made (see also Brami 2015). One potential problem, however, is that evaluating together dates with small and large error margins arising from several generations of radiocarbon dating places too much emphasis on the latter. As already pointed out elsewhere (Brami, Heyd 2011.173), dates of mainland Greece, which were mainly processed in the 1950s -1970s, have two to four times larger standard deviations on average than dates of Western Turkey, making any comparison problematic at best.
The second approach thus incorporated a degree of chronometric hygiene to monitor the quality of the database. A cut-off value of 100 years BP was arbitrarily set for the standard deviation, meaning that 14 C dates with an uncertainty superior or equal to this minimum standard were excluded. Radiocarbon age uncertainty is linked to a variety of factors, not least the resolution of the dating equipment; larger standard deviations may indicate problems with the sample or with the laboratory treatment (Flohr et al. 2015). The problem of 'old wood' effect in charcoal samples was addressed in the following way. First of all, bulk samples, in which carbon of unknown provenience from the sediment is mixed with carbon from the charcoal, were systematically excluded from the audited dataset. Similarly, unidentified charcoal samples which may stem from the inner rings of a tree in which 14 C has started to decay years before the tree was felled or burned were excluded (Zilhão 2001.14181). Finally, long-lived tree charcoal samples from structural timbers such as posts and roof beams which could be reused in successive buildings (Cessford 2001) were flagged out and the corresponding dates discarded.
As a result, short-lived materials such as cereal grains, hazelnut shells, bone/antler made up the essential part of the audited dataset. Bone was treated with caution: bone samples from before the introduction of AMS dating (e.g., four UCLA dates from Argissa) were excluded; likewise, burnt bone and bone apatite (Flohr et al. 2015). This approach is also not without problems. Human bones from coastal regions and river valleys may still have a reservoir effect due to human consumption of marine resources. Seeds, on the other hand, are prone to move across the sediment and, conversely, may be too young. Another consequence is that the dataset on which the second kriging simulation was based was significantly reduced, to 280 dates from 26 sites, leading entire regions such as Greece to be interpolated from only a few known sites. In conclusion, each of the two methods of sampling, selective and non-selective has advantages and limitations, but we argue that, taken together, they provide a valuable snapshot of early agricultural expansion out of Anatolia.
With regard to the input data that was fed into the kriging, it consisted of exact calendar dates (Tab. 1). Since calibrated dates are always expressed as a possible range between two values, not as a specific point in time, a protocol was followed to artificially determine the most statistically probable starting date of each site (Fig. 2). The method of median estimators of phase boundaries was used (Bronk Ramsey 2009a; see Thissen 2010 for a practical application). In short, a Bayesian model was created for each site in which sufficient stratigraphic and contextual information was available for the units sampled (e.g., chronometric phases based on ceramic evidence). Bayesian modelling narrowed down the statistical interval of the dates using prior information about, inter alia, the relationship of the dates, for instance their belonging to the same stratigraphic phase, or their coming 'before' or 'after' one another. In practice, this was done using the boundary function in OxCal 4.2 (Bronk Ramsey 2013). Outliers' dates, showing poor individual agreement (A< 60%) between the observed data and the model, were identified and down-weighted using the outlier analysis approach described by Bronk Ramsey (2009b). A uniform prior probability of 0.05, corresponding to a 1 in 20 probability of each sample being an outlier, was selected (Bronk Ramsey 2009b.5; see also French, Collins 2015.125). Finally, the median was used as the point estimator for the start phase (Thissen 2010).

Kriging interpolation
The dispersal of early farming from Central Anatolia to the Southern Balkans was modelled using the kriging technique of spatial interpolation and the 14 C values derived above. The principle of kriging is that, knowing the value of a set of points in space, it is possible to estimate the value of other points for which data is absent. This is based on the measure of spatial autocorrelation, expressed through a variogram. The variogram is a function describing the degree of spatial dependence of a spatial stochastic process (Wackernagel 2003). Its calculation is based on the distances among the available paired observations. A mathematical model can hence be fitted to the experimental variogram and the coeffi- cients of this model can be used for the estimation through the kriging regression (for more information regarding the statistical process, see Cressie, Wikle 2011). Bocquet-Appel and Demars (2000;Bocquet-Appel et al. 2009) applied this method based on the known distribution of 14 C dates on a uniform grid, in order to estimate the advance of a pioneer front within the context of a colonisation process. This method has some limitations; in particular, it is based on an assumption of spatial homogeneity (Krivoruchko 2012;Pilz, Spock 2007). In other words, this technique appears to be very effective when a subjacent trend is found. Fitting the variogram model to the observed data is a delicate process, which influences the parameters of the regression; if a spatial correlation is not evident, the risk of using an unsuitable variogram model is high.
For the present research, it is not clear from the outset whether or not the data has a linear distribution, so it is hard to find a good predictor for it with ordinary kriging. However, in many cases, the best predictor can be non-linear: empirical Bayesian kriging is a method for predicting non-linear distributions (Krivoruchko 2012). Empirical Bayesian kriging accounts for the error introduced by automatically drawing the variogram trend from a range of individual trends. The new variogram models are estimated on the basis of the previously simulated data; a weight for each variogram is given using the Bayes' rule, showing how likely the observed data could be generated from this variogram. The result of this procedure is the creation of a spectrum of variograms. The predictive density can be calculated by averaging transformed Gaussian distributions (Pilz, Spock 2007).
The variogram for the comprehensive dataset is shown in Figure 3. In order to make the calculation of distances the most accurate possible, the sites are in a metric projection (Universal Transverse Mercator). The values on the x-axis are expressed in metres raised to 10 5 (1 = 100 000m = 100km) and show the distances among the observed points; the y-axis, in turn, shows their semi-variance. The very high variance near the origin indicates a local heterogeneity, added to unavoidable issues related to the 14 C dates themselves (e.g., data quality, dates not belonging to the earliest Neolithic horizon in the region). The low slope of the estimated variograms shows a very low spatial correlation. The variograms for the audited dataset is represented in Figure 4. In this case, the variance at the origin is much lower, and the trend of the simulated variograms shows a higher spatial correlation. Therefore, this dataset appears more appropriate to represent the spread of Neolithic farming. These variograms are inputted in the kriging interpolation model, providing a graphical representation of the possible timing and path of the spread through the use of isochrones, which are boundaries that contain homogenous dates.

Summed probability distributions
In addition to the kriging, the calibrated probability distributions of all 14 C dates were summed in order to gain an insight into regional population fluctuations. This approach rests on the assumption that the density of radiocarbon dates in the dataset is directly proportional to human activity (Steele 2010). In fact, both research and taphonomic biases are likely to affect the shape of the 14 C frequency distribution. To avoid sites being over-represented in the dataset (e.g., Çatalhöyük East alone accounts for over 19.4% of all accepted 14 C dates in the study region), multiple radiocarbon dates for each site were first summed to a single distribution. These distributions were then summed across four target regions (Fig. 5).

Fig. 3. Spectrum of the semivariogram models produced by empirical Bayesian kriging for the comprehensive dataset.
Summed probability distributions in this case may not be used as accurate demographic proxies, given that the number of radiocarbon determinations in each region is below the 500-date minimum threshold quoted in the literature (Williams 2012). This approach, we admit, leaves open many issues; in particular, peaks and troughs in the distribution may not necessarily reflect population expansion and decline, but instead the plateaus and wiggles of the calibration curve (Williams 2012.581). The aim here was to detect major regional discrepancies in the dating, in the order of several hundred years; summed probability distributions provide a valuable medium to show just how well certain periods are represented in terms of 14 C date distribution. They provide an additional control layer, showing not just when farming initially took off, but also how this process was sustained over time, once all the dates are taken into consideration.

Results
The kriging interpolation of the space-time distribution of 14 C dates, whether based on the entire dataset or only a sample thereof, indicates a westward regression of the onset of farming from the Central Anatolian Plateau to the Aegean Basin, followed by a northward shift to inland Thrace and Macedonia. The incremental way in which the isochrones ripple out of Central Anatolia may, we argue, be an artefact of the kriging. Multiple isochrones, at short distances from each other, presumably indicate a standstill or very slow progression. In turn, summed probability distributions of calibrated radiocarbon dates indicate that the advent of farming in Western Anatolia was delayed by up to 2000 calibrated years, supporting the identification of a major chronometric lag between the start of the Neolithic in this region and in Central Anatolia. Figure 6 shows the expansion of the Neolithic, in 250-year isochrones, based on a comprehensive dataset of modelled radiocarbon dates. Compare with Figure 7, which draws on the modelled values of the audited dataset, while sharing the same simulation environment. Both simulations highlight the remarkably early uptake of agricultural production on the Central Anatolian Plateau, which was presumably a major centre of food-plant and animal domestication (Buitenhuis 1997;Asouti, Fairbairn 2002;Martin et al. 2002;Pearson et al. 2007;Arbuckle et al. 2012). Surprisingly, the Pisidian Lake District, which is located at the western end of the Anatolian Plateau, already reflects a much younger tradition. The interpolation shows between two (Fig. 6) and six (Fig. 7) 250-year isochrones between Cappadocia and the Lake District, that is, a little over 200km, in a region which is not characterised by any major topographic boundary. If there was an expansion of the Neolithic towards the west, across the Anatolian Plateau, it was extremely slow-motion, possibly lasting hundreds if not thousands of years. The second kriging simulation, in particular, struggles to interpolate this advance, marked out by not too distant sites showing major discrepancies in corrected start date value, e.g., Asıklı (7934 calBC) and Höyücek (6353 calBC). The kriging produces artificial contour lines to span what is essentially a major lag between two Neolithic regions. In any case, the pattern suggests that agriculture was initially held off in Cappadocia and the Konya Plain, with the 'bond' finally breaking sometime in the 7 th millennium calBC (Düring 2013).

Modelling the advance of the agricultural pioneer front
This above-outlined view is further supported by the subsequent change in direction of the isochrones, from south to north, rather than from east to west, in the Aegean Basin. Here, the two simulations differ significantly. The first kriging simulation based on the non-audited dataset suggests that the Lake District, together with Knossos in Crete, provided a starting point for the initial spread of the Neolithic into Europe. Once the older dates from Hacılar and Bademagacı are excluded from the dataset, due to their poor quality, the Lake District becomes a potential crossroad between a land-way, from the west across the Anatolian Plateau, and a sea-way to the west, spearheaded by slightly older sites like Çukuriçi Höyük and Ulucak. At present, the chronological differences between the Lake District and the Aegean coast of Anatolia are too small to draw firm conclusions about the existence of this second route.
The first kriging simulation highlights a fairly synchronous adoption of agriculture on both sides of the Aegean Basin (Fig. 6). If true, the Aegean Sea probably acted more as a bridge than as a frontier, as also indicated by early dates on the islands of Crete (Knossos), Kythnos (Maroulas) and Gökçeada (Ugurlu). Southern Aegean sites appear to be slightly older than those in the north on average by between c. 500-750 years depending on the simulation, but the distance to cover is much greater, approx. 600km from one end of the Aegean Basin to the other. Once again, differences in the dating are significant but not drastic; they may be explained by other factors, such as a plateau in the calibration curve in the first half of the 7 th millennium calBC, which may influence the simulation (Reingruber, Thissen 2009;Weninger et al. 2014). On the other hand, radiocarbon dates for the Aegean seaboard sites and adjacent regions, like the Thessalian plains, are significantly older than those encountered further inland, particularly in Thrace. Upriver sites in the Struma and Maritsa valleys demonstrate at least one further chronological step in the advance of the Neolithic, with the resulting expansion potentially being driven from west to east rather than from east to west (Lichter 2006).

The Central/Western Anatolian farming frontier
The rapid and incremental manner in which the interpolated isochrones succeed each other across the Central Anatolian Plateau (Fig. 7) lends support to the idea that agriculture was initially contained within this region, spreading internally to multiple sites and communities before radiating outward (Düring 2013). A regional stasis at the onset of the Neolithic in Anatolia can be represented graphically using summed probability plots (Fig. 8A-D). Notice in Figure 8A the calibrated probability distribution of 14 C dates in Central Anatolia during the interval 8500-7000 calBC. Remarkably, this period is almost entirely unaccounted for in Western Anatolia (Fig. 8B),

Fig. 5. Location of the four target regions: A: Central Anatolia (red); B: Western Anatolia (blue); C: Greece (orange); D: Thrace (green). Background map designed by M. Börner.
suggesting that the Neolithic started there between 1500-2000 years later (Brami 2015). The peak in distribution after c. 6500 calBC perhaps marks the initial explosion of the Neolithic in the region. As it is barely noticeable in the other graphs, this peak is unlikely to have been artificially created by the calibration process.
Within the current dataset, there is no indication that Neolithic expansion in Western Anatolia was preceded by a population crash in Central Anatolia. On the contrary, the two distributions run largely concurrent during the interval 6500-6000 calBC. Further west, the question can be raised as to whether the abandonment of sites in Western Anatolia post c. 5800 calBC coincided with a renewed expansion of the Neolithic into Greece and Thrace. The summed probability distribution for the Greek Neolithic is skewed towards a slightly later horizon (Fig.  8C). Dates spanning between c. 7600-7000 calBC are statistical outliers, which can be firmly discounted (Perlès 2001;Brami, Heyd 2011). They show the inherent risk involved in keeping dates with large standard deviations from old excavations generating background noise, as in this case. Thrace represents a further step in time, with the greater part of the distribution presumably falling outside the study period (Fig. 8D).
For reference only, the rate of expansion of the Neolithic for the region under review was measured using the technique described by Ammerman and Cavalli-Sforza (1971). A regression to calculate the rate of expansion, as per the cited article, was performed, using Asıklı as a potential centre of diffusion. The speed implied by the distance-versus-time regression was 0.32 ± 0.11km/year (the range of 0.11 corresponds to the 95% confidence interval), while the time-versus-distance regression returned a much faster diffusion rate, 1.07 ± 0.36km/year (Fig.  9A). The first regression (distance-versus-time) would be preferable if most of the error were due to the dating, while the second (time-versus-distance) if the error were due to the distances. In this case, the distances are exact, so the first regression is of more direct relevance. This approach assumes a linear fit of the regression coefficient. For the present dataset, the correlation coefficient was low, i.e. 0.58 (compare with >0.80 in Pinhasi et al. 2005). This relatively low spatio-temporal correlation is illustrated in figure 9B, where the data distribution appears to be divided into two clusters. Data clusters show the potential lag in Neolithic occupation between Central and Western Anatolia, further undermining the relevance of a linear fit.

The case for an arrhythmic model of Neolithic expansion
If we assume a linear regression from a hypothetical origin in Asıklı, the rate of expansion of the Neolithic within the study region was very low, 0.32km/year on average (Fig. 9). It was much lower, for instance,   (Ammerman, Cavalli-Sforza 1971.681;1984;Pinhasi et al. 2005). Marina Gkiasta et al. have already pointed out that Ammerman and Cavalli-Sforza's average concealed wide regional variations: only 0.7km/year in the Balkans, but a record 5.6km/year in Central Europe (Gkiasta et al. 2003.45;see Ammerman, Cavalli-Sforza 1971. 684). In what follows, we suggest that calculating a mean rate of expansion for the study region is potentially misleading, because it assumes a linear wavedispersal model, which is not consistent with the evidence . From a regional perspective, one indeed observes that linear regression models unduly normalise highly particularised sets of values.
Several arguments can be made in support of an arrhythmic model of Neolithic expansion. First of all, the non-uniform distribution of the isochrones in the two kriging simulations and their change of direction over time (east-to-west then south-to-north) are strong indications that farming did not expand in a linear manner, spreading in fits and starts . Furthermore, the incremental way in which the isochrones ripple out of Central Anatolia in the second simulation (Fig. 7) suggests that farming expansion in this region was extremely slow or halted. A long stasis at the outset of the Neolithic on the central Anatolian Plateau has been represented graphically using summed probability distributions (Fig. 8). Data clusters in the age-distance graphs further de-monstrate the existence of a chronometric lag between Central Anatolia and regions further afield (>400km; Fig. 9).
The results outlined in this paper are consistent with a previous identification of a 2000-year lag in Neolithic occupation between the central Anatolian Plateau and the Aegean Basin (Brami 2015). Farmers appear to have been initially held off in this region. On account of the summed probability plots, there is no indication that a 'bust' preceded the 'boom', as in other regions of Europe (Shennan et al. 2013). No regional population collapse can be detected in Central Anatolia before c. 6000 calBC (Fig. 8A). On the face of the evidence presented, the idea of a farming frontier crystallising as a result of either a loss of momentum in the Neolithic core or an encounter of resistance in Western Anatolia appears more likely. The 'bond' was finally breached c. 6500 calBC, with a subsequent explosion of sites recorded throughout Western Anatolia (Düring 2013).

Limitations of the study
Kriging is arguably a powerful technique to interpolate the spread of early farming across Eurasia (Bocquet-Appel et al. 2009). One issue that this paper has sought to address is the assumed linearity of ordinary kriging, which makes the computation of non-linear expansion behaviour, such as an arrhythmic spread in fits and starts, problematic. Where sites on either side of a 'frontier' display widely different values, ordinary kriging breaks down the gap be-  tween them into a series of isochrones, essentially imposing linearity where there is none. The method of kriging which was used here, empirical Bayesian kriging addresses this issue by adjusting the simulation at each of the input data locations (Krivoruchko 2012). Although the results obtained with this method indicate an improvement in kriging data with non-stationary covariance structure, the second interpolated map (Fig. 7) still displays an incremental pattern of expansion out of the Central Anatolian Plateau ('ripple' effect). One potential issue with this simulation lies in the number of plotted sites, which at 26 is not high enough to generate an accurate isochrone map. The second kriging simulation possibly lacks in resolution what it makes up for in data quality.
Another limitation of the kriging method of interpolation as it has been pursued here is that it operates in a spatially neutral environment, where every section of the map is given equal weighting regardless of its geographic context, i.e. valley bottom, mountain top, sea, etc. This is consistent with previous applications of kriging for modelling the expansion of the Neolithic in Eurasia (Bocquet-Appel et al. 2009).
One way forward would be to use the 'best patch' variable (e.g., Bocquet-Appel et al. 2014.63-64). This amounts to grading land according to their agricultural potential. Approaches based on the spatiotemporal distribution of 14 C dates are helpful to describe a geographic spread, less so to analyse or explain it. The models presented in this paper do not take into account a multitude of variables which may have influenced early farmers. A different approach, which estimates climatic variables and their effect on the landscape as well as the socio-economic systems and demographic structure, is agent-based modelling. This holistic approach, which brings in data from different disciplines (economy, anthropology, ethnography, paleo-climatology), has been recently introduced in archaeology, allowing one to test scenarios that could not be inferred from purely archaeological observations (Axtell et al. 2002;Kohler et al. 2007;Janssen 2009;Bocquet-Appel et al. 2015).

Conclusion
This article has established, through a suite of geostatistical and graphical simulations, that the advance of the Neolithic from the Anatolian heartland to Southeast Europe involved at least two distinct stages. Farming was initially held off on the central Anatolian Plateau. Up to 2000 calibrated years were necessary to bridge the chronometric lag between Central and Western Anatolia (Brami 2015). Once early farming finally spread into the Southwest Anatolian Lakes Region and the Aegean Basin, shortly before c. 6500 calBC, it rapidly made its way north, reaching Eastern Thrace c. 6000 calBC. The pattern of spread described in this paper is consistent with an arrhythmic model of diffusion, involving major standstills (or 'arrhythmic phases')i.e. the Central/ Western Anatolian farming frontier -punctuated by rapid and/or regular advances in the Aegean Basin and the Southern Balkans (Guilaine 2000.268-270).
Moreover, this paper has demonstrated that linear regression models, such as the 'wave of advance' (Ammerman, Cavalli Sforza 1971;1984), virtually conceal strong regional variations in the data by normalising them. While these approaches may be useful on the scale of Eurasia to describe the overall pattern and direction of spread, moving one scale down, the reader can see that they fail to reflect the fits and starts of the process, in this case the crystallization of boundary or frontier zones, which preceded the ultimate explosion of farming communities c. 6500 calBC. The use of Bayesian clustering alongside the kriging helped to further sharpen the resolution of the model. Two kriging simulations were presented, one based on virtually all calibrat-ed 14 C dates, the other on a strictly audited sample. Together they provided a valuable picture of the Neolithic expansion out of Anatolia, as evidenced by the time-space distribution of 14 C dates. Finally, summed probability plots were used to show regional population fluctuations past the initial expansion of the Neolithic.

M. Brami is the recipient of a post-doctoral fellowship of the National Research Fund, Luxembourg (AFR project no 9198128). A. Zanotti is a BEAN researcher (Bridging the European and Anatolian Neolithic).
Thanks are owed to Barbara Horejs, Bernhard Weninger and to members of the Anatolian Aegean Prehistoric Phenomena (AAPP) research group at OREA for helpful information and comments. The responsibility for any shortcomings remains with the authors.