The origins of agriculture in Iberia > a computational model

Salvador Pardo Gordó1, Joan Bernabeu Aubán 2, Oreto García Puchol 3, C. Michael Barton 4 and Sean M. Bergin 5 1 Departament de Prehistòria i Arqueologia, Universitat de València, ES salvador.pardo@uv.es 2 Departament de Prehistòria i Arqueologia, Universitat de València, ES jbauban@uv.es 3 Investigadora Ramón y Cajal\Departament de Prehistòria i Arqueologia, Universitat de València, ES oreto.garcia@uv.es 4 Center for Social Dynamics and Complexity, Arizona State University, USA michael.barton@asu.edu 5 School of Human Evolution and Social Change, Arizona State University, USA sean.bergin@asu.edu

in research on the spread of agriculture.In this context, computer simulation has become one of the techniques most frequently used to explore the space/time of Neolithic dispersal and its subsequent evolution.
The introduction of computer applications in archaeological research can be dated roughly to the 1950s.The first work focusing on simulation per se was Doran's short essay on cybernetics and its application as a useful tool for generating explanations of the archaeological record (Doran 1970.296-298).Subsequently, computer simulation applied to the problem of the dispersal of the Neolithic can been found throughout the archaeological literature for over 40 years.The first and most influential work was framed by Albert J. Ammerman and Luigi L. Cavalli-Sforza (1971;1973;1979;1984), which was based on an adaptation of Fisher's reaction-diffusion model applied to the spread of agricultural groups driven by constant population pressure, so-called logistic growth.They evaluated this model for the diffusion of agriculture across different areas of western Eurasia (1984.134-135)by comparing the timing of the initial arrival of agriculture predicted by their model with then-available radiocarbon dates from the archaeological record.They concluded that the predictions of their model and the archaeological information strongly correlated (R ~0.8).They also suggested a southeast-northwest gradient for the spread of agriculture across Europe, validating the theory of a Near Eastern origin for the Neolithic as promulgated by Grahame Clark (1965).Although we are discussing the Neolithic expansion in Europe here, other simulation work has focused on the spread of rice in Asia (Silva et al. 2015), and the expansion of Paleolithic populations (Fort et al. 2004) or languages, such as Bantu (e.g., Grollemund et al. 2015;Russell et al. 2014).
In the past 15 years, the availability of inexpensive, high-speed computer processing and a greatly expanded radiocarbon database has led to a number of studies revisiting the empirical comparisons and demic diffusion models of Ammerman and Cavalli-Sforza, using different approaches such as timedelay, the role of waterways, effects of boundaries and cultural practices (e.g., Ackland et al. 2007; Davison et al. 2006;Fort et al. 2012;Fort, Méndez 1999).In other research, we conducted a detailed review of some of the most notable such work (e.g., Bocquet-Appel et al. 2009;Davison et al. 2009;Gkiasta et al. 2003;Pinhasi et al. 2005), concluding that new radiometric information from the Iberian peninsula has not yet been fully utilised in computer models for Neolithic dispersal at continental scales (Pardo Gordó et al. in press).This large body of new radiocarbon dates only has been used in local spreading models (Bernabeu et al. 2015;Isern et al. 2014).
Since the 2000s we are now in a position to highlight the growing interest in examining different theoretical frameworks by means of archaeological simulation, and the corresponding increase in the number of papers focused on modelling work (Costopoulos 2010; Lake 2014).Computational modelling has become a more common and sophisticated tool in the archaeological analytic toolbox (Barton 2013a;2013b), although the use of computers to support social theory more generally is hardly actu-

Fig. 1. Map of the Iberian Peninsula with Early Neolithic sites with radiocarbon dates used for the model evaluation.
ally a new concept (Hägerstrand 1965).In this paper, we investigate the spread of agriculture in Iberia using by means of simulation methods, and compare results with the preliminary models for the spread of the Late Mesolithic, the so-called Geometric Mesolithic.We focus on the Iberian Peninsula because it is a particularly good region in which to study the process of agricultural dispersal.It has evidence of populations of foragers during the final Mesolithic, post quem 6000 BC (Bernabeu et al. 2014).It is situated at the western extreme of the Mediterranean Basin and serves as a bridge between Africa and Europe.For these reasons, Iberia can be considered a sub-continent where it is possible to examine a number of processes related to the Neolithic transition.For example, this area is the best place to evaluate the possibility of duel expansion routes (South-eastern France and Northern Africa) of the first groups of farmers.This has become a topic of interest recently, although there are different views on its impact on the process of Neolithic expansion (see Cortés Sánches et al. 2012;García Borja et al. 2014;Zilhão 2014 for references).

Computational model
We use computer simulation models, more specifically in Agent-based Model (ABM), to investigate the spread of agriculture in Iberia.This methodological approach is one of the most active applications of simulation in archaeology (Lake 2015) despite its lack of use in studies of the spread of farming (Parisi et al. 2008).Briefly, ABM is a kind of computational model with agents that are discrete and autonomous entities that differ from others in space and time, and usually interact with others or with their environment locally (Bonabeau 2002;Railsback, Grimm 2012).
Our spread model (Bergin et al. 2015) was implemented the Netlogo modeling platform (Wilensky 1999) because it allows us to import and use geo-referenced datasets within the modelling environment, including radiocarbon dates and other kinds of information (in our case, ecological).For this reason, our model takes the form of a spatially explicit cellular automaton in a gridded landscape in which agriculture can spread on the basis of rules of dispersal.Our approach is based on "modelling as experiment" (Bankes et al. 2002) as this allows us to use computational model environments to explore the effects of different variables and compare hypotheses to existing datasets (Grimm et al. 2005).

Virtual world
Currently, the emphasis on the importance of environmental conditions is a triggering factor for the dispersal of Neolithic groups (Gronenborn 2009;2010).Although it is widely recognised that ecological contexts are more or less suitable for early Neolithic agriculture, this has not been considered explicitly -with a few exceptions -in the modeling work (e.g., Ackland et al. 2007;Banks et al. 2013).
We classified landscape cells based on their suitability for cereal agriculture, using a combination of terrain and climate parameters 1 1 (Bevan, Conolly 2004;López Bellido 1991).We focused on wheat, because it has the most stringent climatic requirements of the different species of early Eurasian cereals.Maps for minimum temperatures for March, maximum temperatures for the spring months of March through May, and total precipitation for spring months were combined to create an index map of suitability for cereal agriculture; these are summarised in Table 1.A combined ecological suitability index was created by summing the three climate index maps and slope index map.The resulting map was scaled to a 5 x 5km resolution and uploaded to NetLogo.Each patch in the models then has a suitability index value based on a combination of the variables described above.

Spread movement, demographic effects and starting points for agriculture dispersal
The three modes of Neolithic dispersal tested in our model are neighbourhood, leapfrog and the Ideal Despotic Distribution (IDD) model (Fig. 2).The first corresponds to the classical wave-of-advance movement promulgated by Fischer (1937) and applied to population expansion by many researchers (see Steele 2009 for references).The model is straightforward: agriculture spreads from one cell to neighbouring cells that lack agriculture as long as they are suitable for it (i.e. have a sufficiently high ecological suitability index value).
The second corresponds to the leapfrog model described by Tjeerd Van Andel and Curtis Runnels (1995).This algorithm simulates the dispersal of agriculture from any cell that has agriculture to another randomly selected cell within a given distance (specified by the user) which does not yet have agriculture and that is suitable.This punctuated spread is also the kind of movement proposed in the maritime pioneers models (e.g., Dawson 2011;Zilhão 2001).Two related types are "neighbourhood with no ecological constraints" and "leapfrog with no ecological constraints".These work like the constrained versions already described, but without taking into account the suitability of cells for agriculture.
The third process is the IDD model from Human Behavior Ecology (Kennett, Winterhalder 2006;Smith 1992;Smith, Winterhalder 2003), it was implemented as a follow-up on suggestions by Stephen Shennan (2008) and Sarah B. McClure et al. (2006) about the potential impacts of socially mediated access to re-sources during the Neolithic.In this case, agriculture spreads to the neighbouring cells with the highest suitability values, but this suitability is affected by the number of farmers already occupying the cell.
That is, values decline whenever agriculture 'spreads' to a cell in which it is already present, and agriculture will spread only to neighbouring cells with the highest suitability values.
Finally, in this model, we explored 17 different potential starting points for the spread of the Neolithic

Fig. 2. Examples of spread models in action. A: shows wave-of-advance dispersal; B: shows the IDD spread algorithm; C: shows leapfrog dispersal with the maximum leap distance set to 5 cells. On the maps, an 'X' marks the starting point for the spread; yellow dots show the locations of Neolithic sites. The colours indicate the relative time of arrival of agriculture: the darkest red is the oldest arrival time, and lightest pink the most recent arrival time. Underlying green shades show the ecological suitability of cereal farming.
across Iberia.We chose the mouths of various rivers or areas near of them (e.g., Málaga and Gibraltar) around the perimeter of the Iberian Peninsula, with one of them in the centre as a null case (Madrid).

Previous results
To estimate a chronological range sufficient to encompass the spread of agriculture over much of the Peninsula, we first identified the oldest acceptable unquestionable date for the use of domesticates: a date of 7569±48 calBP (all dates used here are expressed as calibrated years BP.)We then extended this range up to 6000 calBP to encompass the earliest evidence for agro-pastoral systems across the Peninsula.This range permits us to cover a total time span of between 7800-6000 calBP, with the last 500 years for sites located only in northern Spain.For any region in the Iberian Peninsula, we selected sites representing the earliest dated evidence for domestic plants and/or animals.The radiocarbon dataset (Bernabeu et al. 2015.Tab. 2 SI) includes only dates clearly associated with archaeological remains of domestic taxa (plants or animals).In total, we have 134 radiocarbon dates associated with 115 archaeological sites.Their distribution can be seen in Figure 1.In total, 53 refer to long-lived taxa, 39 to shortlived taxa and 42 to domestic taxa (Fig. 3).We grouped this radiocarbon information into four subsets (the mean radiocarbon age is used in all groups): ❶ Best: includes a mix of dates made on domestic taxa where available, non-domestic short-lived taxa when directly dated domestic taxa are not available, and non-domestic long-lived taxa when this is the only kind of radiocarbon sample available.Our first work (Bernabeu et al. 2015) focused on exploring the radiometric dating sample, points of origin for the Iberian Neolithic and exploration of parameters such as movement, distance, ecology and occupation costs.In the first experiment, we evaluated archaeological samples and initial expansion points, keeping the values of movement, distance and cost of occupation fixed (Bernabeu et al. 2015.Tab. 1).The results show that the samples used influence the results, and the best starting points are systematically located in eastern Spain, confirming the Mediterranean origin of the Neolithic.In the second experiment, we evaluated whether the fit between the model and the empirical data improves with multiple origin points instead of a single origin point.This experiment allowed us to test a possible double entry route for the Iberian Neolithic.The results of this experiment allowed us to discard the idea that simply increasing the number of origin points increases the correlation results.
We concluded that 9 of the 10 strongest correlations are associated with a dual entry route of the Neolithic into Iberia (one of them located in the northeast and the other in the southeast) and a complex, multispreading process.
Finally, using the best correlations of the previous experiments, we explore movement, distance, ecology threshold and the costs of existing occupation by farming groups.We observed the best correlations are associated with leapfrog dispersal, with a distance between 25-50km, medium-high impacts of prior agricultural occupation (demographic aspects) and a preference for places with high potential cereal productivity (ecological threshold between 5 and 6).This allowed us to conclude that the expansion of Neolithic into Iberia can be characterised by pioneer colonisation, whereby farmers travelled relatively long distances looking for places with no or few people already farming, and an attractive environment for wheat.
Finally, in other work (Pardo Gordó et al. in press), we explored in more detail the radiocarbon data and its influence on our model results with several experiments.The first compared different groups (above) from the radiocarbon dataset, with a single origin point, and more specifically the best and oldest sub-sets.We observed that that 15 of the 20 strongest correlations are associated with the best sub-set, suggesting that different selections of the radiocarbon information can produce quite different results.Next, we compared the best sub-set with short-lived dates.Again, we looked at the 20 strongest correlations, with unexpected results.The more 'reliable' short-lived radiocarbon dataset generated correlation coefficients considerably worse than the larger, mixed best dates set.Why?We conducted a sub-experiment to test whether dated shell that had potentially been affected by the reservoir effect (Ascough et al. 2005;Soares, Dias 2006) could have had an impact on the results.We again selected one starting point (the Segura River, eastern Iberia) for each of the 5 configurations and removed those dates for shells in the short-lived data set.Removing shell dates from this sub-set significantly improved its match with model results.It is worth remembering that the use of samples made on shells can be problematic when used to evaluate model results if the reservoir effect is not taken into consideration.
In the last experiment, we compared the short-lived dates with the smaller group of dates from domestic taxa.Of the 25 best correlations, better Pearson correlations coefficient were produced from the more reliable dates of domestic taxa only dates than the larger short-lived dataset, even without dates for shell.
In short, our previous work suggests that the quality of the radiocarbon information used needs to be considered carefully when using a body of dates to evaluate the results of computational modelling of the spread of farming (empirical evidence for this new economy).The importance of using careful and rigorous criteria for the selection of radiocarbon dates noted by other archaeologists (e.g., Bernabeu 2006;Zilhão 2001;1993;2011;Bernabeu et al. 2001;Bernabeu, Martí 2014;Rojo et al. 2008) is firmly reflected in the results of our modelling experiments.Nevertheless, the poor results obtained from samples made on short-lived taxa associated with domestication economies were surprising.

Auditing radiocarbon problems, new modelling results
As we observed in the section above, the best correlations obtained from previous experiments made on remains of domestic and dates on short taxa (including domestic and non-domestic plants and animals), generated Pearson correlation coefficients considerably worse than other subsets including the oldest and the best.We suggested that these poor correlations could relate to the reservoir effect (on shells and bones).Consequently, we need to calculate the reservoir effect and its impact on spatio-temporal variations (for details see Ascough et al. 2005).
As we pointed out (Bernabeu et al. 2014), these problems are especially visible in Portugal, where a significant number of dates derive from shells and human bones.
Also, as recently pointed out by Rachel Wood (2015) and Karl-Göran Sjögren (2011), problems linked with the sampling criteria can also affect different treatment procedures in the laboratory.At the same time, the ratio of nitrogen to carbon in bone collagen has been proposed as a good indicator for testing the quality of radiocarbon results (Van Klinken 1999).Unfortunately, the details of the N/C ratio are not usually available for the published radiocarbon dates, adding uncertainty about the possible importance that this kind of problem in radiocarbon assays of bones.Finally, Haidé Martins and colleagues (2015) demonstrated that distinguishing some domestic taxa in animal bones (especially Ovis sp. in the Iberian Peninsula) can be difficult, with consequences for dating the beginning of farming.Bearing in mind the potential effect in the radiocarbon outputs, we designed a new experiment that considers only charred samples such as seeds, fruits and char-coals identified as short taxa (shrubs) and we add domestic bones only when the N/C ratio is known and adequate for dating.A total of 34 radiocarbon dates meet these criteria and were used for the experiments reported here (Tab.2). Figure 4 shows a comparison of our previous results obtained with domestic taxa (for details see Bernabeu et al. 2015. Tab. 1) and the results using the same model parameters obtained using the new audited radiocarbon data set.As shown in the graph, the correlation obtained increases significantly.
To further illustrate this point, if we look at the results associated with the point of origin set to the Rio Segura and using the wave-of-advance spread algorithm with ecology considered, the use of domestic taxa shows only a value of R = -0.39,while the use of a database with the filtered information increases its correlation to R = -0.50.
In sum, these results suggest again that the radiocarbon samples used have significant effects on the correlations obtained, and consequently on the evaluation of different model scenarios.If we want to be sure about the evaluation of our models (including mathematical, agent-based or cellular automata) to analyse Neolithic dispersals (and, of course, other similar phenomena) using radiocarbon dates, then we need to carefully audit the samples, a task on which we are working now in order to reexamine our previous conclusions (Bernabeu et al. 2015;Pardo Gordó et al. in press).

Geometric spread as a null hypothesis
Mesolithic bladelet technology, including trapezoidal forms appeared in the 9 th millennium calBP as a European phenomenon which included the appearance of new techniques and tools in lithic industries.A millennium later, agriculture expanded around Western Europe.The Mesolithic dispersal has been considered by several authors, such as Clark (1958), who compares this expansion with the posterior Neolithic advance.Despite an interest in exploring the mechanisms behind this dispersal (demic versus cultural), only a few works have highlighted this potential line of research, without developing it further (Binder et al. 2012).Instead, most authors focus on the geographical origin of the Mesolithic, arguing over the different potential starting points (Biagi, Kiosak 2010;Binder et al. 2012;Marchand, Perrin 2015).Although there is broad spatial variability in Mesolithic technology across Europe, it is generally thought to indicate a major shift in blade technology and the production of compound arrowheads (geometric tools).This involves knapping techniques to obtain regular blades and bladelets using indirect percussion or pressure as a distinctive characteristic in order to make regular blades for geometric forms (trapezes) with symmetric or asymmetric shapes (Binder et al. 2012).Other tools, such as notched blades, are also common, and were probably used for processing plant materials (Gassin et al. 2013).
In the Western Mediterranean, this cultural complex is known as the Tardenosien tradition, or referred to as the Late Mesolithic.This encompasses the re- gionally named industries of the Castelnovien Complex (or Second Mesolithic in France and Italy), the Upper Capsian in North Africa (Rahmani 2003) and Geometric Mesolithic in Mediterranean Iberia and Portugal (Fortea 1973;Utrilla, Montes 2009).With some regional particularities, this Mesolithic phenomenon has been considered to have across spread Europe in some kind or diffusion process (Kozłowski 2009).
Building on our prior work, we selected radiocarbon dates corresponding to the first Geometric Mesolithic in order to compare some parameters related to Mesolithic and Neolithic dispersals.Current information shows that the Late Mesolithic is well documented in eastern Iberia and the Ebro valley (Mediterranean region), and central and southern Portugal (Atlantic coast).While several authors consider some settlements in the Cantabrian region as Mesolithic with geometrics (Arias, Fano 2009), these settlements did do not include all of the technological elements of the well-defined Late Mesolithic of the Castelnovien tradition, so they were eliminated from our database for this preliminary assessment.Other areas (northeastern Iberia in Catalonia and the inner territories of the Meseta) lack archaeological data on this period.
We compiled a total of 21 dates associated with Mesolithic contexts, considering only audited short-lived samples as described above (Tab.2).The criteria followed the protocols used in our previous work (Bernabeu et al. 2015), considering the most ancient date for each site provided by short-lived samples and comparing them with the modeling results.A particularity in relation to the nature of the samples affects Portuguese Mesolithic contexts, where human skeletons constitute the main material dated.For this, we used the radiocarbon dates compiled by António Faustino Carvalho (2010).
In this experiment, we compare different starting points for the spread of the geometric tools around the perimeter of Iberia and evaluate the modelling results against radiocarbon dates made on shortlived taxa.The parameters for this experiment were set as follows: threshold for ecological suitability (i.e. for wheat cultivation) 0 and 3, costs of prior occupation 5% and leapfrog radius distance of 5 cells (25km).As we can see in Figure 5 that the best correlation between the model result and dated Late Mesolithic sites occurs when the ecological threshold is limited to 0 with R = -0.32 in the best case.
Regarding the best correlations (those that have negative values), we note several results.First, most of the points of origin with negative correlations (except Bilbao) are located on the Mediterranean coast of the Iberian Peninsula.These results parallel the proposed expansion of the Mesolithic complex throughout Europe (e.g., Clark 1958).The best fitting spread algorithm in all cases is the wave-ofadvance (spreading to neighbouring cells only), and when ecological suitability is not considered.
However, are there any similarities between these results and those related to the first groups of farmers? Figure 6 shows the comparison between the Mesolithic and Neolithic (using only dates from domestic taxa).The graph shows that the correlations associated with the Neolithic are higher than those for the Mesolithic, and that the best Neolithic correlations (R > -0.3) are associated with scenarios where ecological suitability is taken into consideration.These results do not seem unreasonable, because the base map used was drafted following ecological parameters for cultivating wheat (see section 2.1), which should not be relevant to Mesolithic foragers.Nevertheless, this first attempt to model the

Concluding remarks
In this paper, we illustrate the potential of bottomup modelling for investigating the dispersal of agropastoral economies and life ways in Europe, focusing on the Iberian Peninsula as a case study.Additionally, we use computational modelling approach as a method of formalising and testing multiple (and complex) hypotheses about local-scale decision rules, rather than as a means of quantitatively characterising agricultural dispersals at the continental scale (so-called top-down models).Agent-based models and mathematical models are complementary approaches to formalising hypotheses about the dynamics of human societies.Top-down modelling allows us to describe general trends and to aggregate behaviour(s) in societies at large scales and over extended periods.On the other hand, bottom-up modelling is particularly well suited to understanding individual behaviour and its interactions with others and its environment, which generated the general trends observed.We believe that the formalisation in both kinds of modelling approaches is an essential step for the ability to systematically compare and test hypotheses about spatiotemporal dynamics of past human societies against a poor, fragmentary and incomplete archaeological record.In short, this paper is a good example of methods useful for understanding a complex problem (the Neolithic spread) with a promising new approach (agent-based models).
Finally, this work demonstrates the importance of carefully auditing the radiocarbon information used to evaluate quantitative models of Neolithic (and others) dispersals.This is essential if we aim to test the reliability of models of human dynamics against the empirical record.

Fig. 3 .
Fig. 3. Bar chart with the number of radiocarbon dates made on long-lived taxa, short-lived taxa and direct taxa.See the online version to identify the colours of each category.

Fig. 4 .
Fig. 4. Correlation coefficients for the results of the Neolithic audited and not audited for individual starting points for agricultural dispersals.The colours indicate the different strategies employed by agents.Positive correlations and models are excluded.

Fig. 5 .
Fig. 5. Correlation coefficients for the results of the Late Mesolithic for individual starting points for agricultural dispersals.The colours indicate the different strategies employed by agents.A red circle indicates negative correlations.

Fig. 6 .
Fig. 6.Correlation coefficients for the results of Late Mesolithic and Neolithic results (only dates on domestic taxa used for comparison) for individual starting points for agricultural dispersals.The colours indicate the different strategies employed by agents.A red line indicates negative correlations > -0.3.

Tab. 1. Environmental parameters used to calculate Ecolo- gical Suitability Index.
In other words, this is the best radiocarbon date for each site.❷ Oldest: the oldest date for each site regardless of the kind of sample.
(Bernabeu et al. 2015;ous results(Bernabeu et al. 2015; Pardo Gordó et al. in  press), we first describe how we compare the model results with the archaeological information.This involves establishing a tem-