Gregarious behavior, human colonization and social differentiation An agent-based model

: Studies of colonization processes in past human societies often use a standard population model in which population is represented as a single quantity. Real populations in these processes, however, are structured with internal classes or stages, and classes are sometimes created based on social differentiation. In this present work, information about the colonization of Old Providence Island was used to create an agent-based model of the colonization process in a heterogeneous environment for a population with social differentiation. Agents were socially divided into two classes and modeled with dissimilar spatial clustering preferences. The modelandsimulationsassessedtheimportanceofgregariousbehaviorforcolonizationprocessesconductedinheterogeneousenvironmentsbysocially-differentiatedpopulations.Resultssuggestthatintheseconditions,thecolonizationprocessstartswithanagentclusterinthelargestandmostsuitablearea.Thespatialdistribu-tionofagentsmaintainedatendencytowardrandomnessassimulationtimeincreased,evenwhengregarious-nessvaluesincreased.Themostconspicuouseffectsinagentclusteringwereproducedbytheinitialconditionsandbehavioraladaptationsthatincreasedtheagentcapacitytoaccessmoreresourcesandthelikelihoodofgregariousness.Theapproachpresentedherecouldbeusedtoanalyzepasthumancolonizationeventsorsupportlong-termconceptualdesignoffuturehumancolonizationprocesseswithsmallsocialformationsintounfamiliaranduninhabitedenvironments.


Introduction
. Human colonization processes are a common research topic in disciplines such as archaeology and population genetics that study social human evolution in our past. Studies of colonization processes in past human societies o en use a standard population model in which population is represented as a single class of agents that o en exhibits highly mobility without strong spatial constraints (Barton et al. ; Jochim ; Winterhalder et al. ). However, real populations during colonization processes are o en structured with internal classes or stages based on age or other natural phases (Briggs et al. ; Hey & Machado ), and classes are sometimes created based on social di erentiation (Rockman & Steele ; Stein ). In this paper we model a colonization process with a socially-di erentiated population in order to provide insights about how the sum of individual preferences and the internal diversity of populations shape human colonization processes. This approach could be used to analyze past human colonization events or support long-term conceptual design of future human colonization processes with small social formations into unfamiliar and uninhabited environments (Billingham et al. ; Smith & Davies ). .
In the work herein, a case study of the European colonization of the Americas was used as input data to model how structured human populations colonize a broad range of environments. A er the arrival of Columbus, Europeans controlled communities and territories of Amerindian societies but also occupied new and unfamiliar landscapes through processes of exploration, migration, and colonization. These processes were undertaken by hierarchically-organized, slave-based, sedentary, and small populations, that stemmed from demographically larger societies with original settlements located in Old World environments. Europeans and enslaved Africans found in the Americas a landscape full of unfamiliar and heterogeneous environments, in which their social conditions and available communication and transportation technologies inhibited immediate interactions with their origins. .
The colonization of Old Providence Island represents one of these cases. Old Providence Island is about km East of the Nicaraguan coast, and km Southwest of Jamaica (Newton ; Ordahl ; Parsons ) ( Figure ). The events that occurred during colonial times in Old Providence Island provide us an opportunity to explore fundamental questions about how humans interact in heterogeneous environments. In this paper, a colonization process was modeled based on information about the colonization of Old Providence Island using agent-based model techniques. Agents were socially di erentiated and modeled with dissimilar spatial clustering preferences. The purpose was to evaluate the importance of gregarious behavior when structured human populations with di erent preferences about group organization colonize a new heterogeneous environment.

Figure : General Location Old Providence Island (a er Esri
).
. The model and simulations were not designed to accurately account for the demographic changes and past land use of Old Providence Island. Instead, we aimed to build what Gilbert calls a middle-range model (Gilbert ). We were interested in exploring how the spatial distribution of agents changes with di erent degrees of independence or interdependence when involved in a colonization endeavor. Independence and interdependence are expected to have an e ect on spatial distribution since they o en require some degree of face-toface interaction (Drennan & Peterson ). Therefore, minimizing or maximizing spatial distance is the most common strategy that humans use to manage costs and benefits of interaction (Drennan ). Agent-based models can provide a digital representation of populations with internal di erences in preference for frequent peer interaction. .
By means of an agent-based model, we simulated a continuum of scenarios between two extremes: one representing a society with a large proportion of independent agents, and the other a society with agents prone to being more interdependent. Agent gregariousness was taken as proxy for interdependence; agents prone to be interdependent were defined as more gregarious while more independent agents were assumed as less gregarious. Scenarios were run in a digital representation of the island's environment. .
The environmental variability and occupation history of Old Providence Island were used as input data to model colonization dynamics. The digital environment was structured with a land classification of the island that used biotic and abiotic data with the aim of identifying areas that were more suitable for the first settlers (Fajardo & Rodríguez ). Social di erentiation was included in the model using the behavioral component of the slavery process on the island. The earliest permanent settlements were built in the th century when English puritans were expanding in the Caribbean. The absence of Amerindians on the island since the beginning of the colonization process led Europeans to bring Africans as slaves. Regardless of the changes in territory control between Englishmen and Spaniards, the presence of Africans was a constant factor in the social organization of the island (Parsons , p. -). An important aspect of the behavior of Africans on the island was that they marooned and sometimes even managed to escape from the island (Newton , -). The voluntary preference to become marooned, namely withdrawing or reducing the possibility of interaction with other individuals by physical migration, was considered as a relevant feature of gregariousness and it was included in the model. This preference was also an important part of the historical context of the island. It highlights how dramatically segregated individuals reacted and decided to act based on their social conditions. Cultural studies underscore how African presence is now in every aspect of Caribbean societies and culture (Hall ). This characteristic suggests that despite their disadvantageous position, enslaved Africans played an active part in the construction of societies created in the Americas. The colonization process of Old Providence was managed by Europeans, but the resulting communities and patterns of social interaction were also shaped by African behavior. For these reasons, the colonization process is referred here as managed by Europeans although the resultant communities are understood as an Euro-African product.

Agent-based modeling in archaeology
. Computer modeling and simulation in archaeology date back to the early s and have explored several aspects of a wide range of societies. Readers can refer to Lake ( ), BarcelÃş et al. ( ) and Grosman ( ) for how computer technologies changed the way archaeologists model, document, visualize, and analyze archaeological data. Recent works also discuss epistemological issues of agent-based computational modeling in archaeology (Saqalli & Van der Linden ; Wurzer et al. ). Lake ( ) highlights that some agent-based approaches with an archaeological scope can be classified as attempts to analyze long-term socio-ecological change as a complex adaptive system, and this is the approach presented here. This is also a topic actively studied in recent agent-based simulations (Lake ; Saqalli & Van der Linden ). .
Complex adaptive system approaches in archaeology that work with agent-based simulations study di erent aspects of agent interactions and settlement systems of pre-industrial societies. These simulations usually develop models focused on understanding one particular long-term sequence of social change. For example, Kohler and collaborators have developed their model as a programmatic research about pre-European populations in the Mesa Verde region of Southwest Colorado. They started from a simple model that applied fixed agent decision-making rules to settlement locations (Kohler & Carr ). Today, their model has grown in complexity and contains algorithms for a decision-making process that includes cultural preferences, kin relations, and a changing environment (Kohler et al. , ). Barton et al. ( ) have been studying agro-pastoralist eco-dynamics in the Mediterranean in order to answer questions about the sustainability of socio-ecological systems. Small bronze age communities from the Middle East have been studied with agent-based models and archaeological datasets that demonstrate resource gains from some household-like agents at the expense of others, producing more economically di erentiated communities through time at the onset of urbanization (Wilkinson et al. ). Other recent evolutionary topics addressed by agent-based simulations have been the process of fission and fusion of settlements in pre-industrial societies (Crema ; Gri in & Stanish ), longterm anthropogenic e ects on landscapes (Barton et al. ), and social cooperation within high risk environments (Clark & Crabtree ).

.
The examples above share an interest in long-term societal change; inclination toward environmental realism as a pivotal part of their model; explicit concern about human-environment interaction; model testing against archaeological and/or historic data; and increasing complexity built over an initial simple model. Although the model presented here shares similar interests in the analysis of long-term human complex systems, it follows a complementary approach. Environmental heterogeneity is included in the model without dynamics. Instead, the model focuses on parameters that represent di erent preferences of social interaction. Also, the analyses do not aim to explain a particular long-term sequence of social change but to explore how diversity in individual preferences about interaction could a ect human clustering in a colonization process.
Aim of the study .
The agent-based simulation was focused on modeling the variability between two simple âĂĲwhat ifâĂİ scenarios. On one hand, a scenario where agents prefer to live in gregariousness near other similar agents. This was taken as a proxy for a population with preference for frequent interaction and dependence on each other (interdependent scenario). On the other hand, a scenario where the majority of agents preferred to live far from similar agents. This scenario represented a population with preference for self-sustainability and low daily interaction (independent scenario). The rationale behind this approach was to simulate the heterogeneity of spatial patterns, based on di erent preferences that agents in a given society could have in their attitudes towards living in proximity to others. . When Amerindian-European-African dynamics in the Caribbean started, Old Providence Island functioned as an uninhabited territory not only for European-African communities but also for Amerindian populations.
By the year CE, Englishmen, women with some children and slaves were living on Old Providence Island. This population was divided in two types of settlements: one small nucleated settlement and several dispersed farmsteads. The small nucleated settlement was named New Westminster and it was located in the Bowden area (Figure , northwest of the island). New Westminster only had houses, a church and the governor's house (Newton , pp. -), while the rest (most) of the population lived outside of New Westminster in dispersed farmsteads. The enslaved African population was introduced in CE. Some Africans managed to escape from landlords and established shelters in forest areas, especially in the southeast part of the island. These shelters were burned down by Europeans at least two or three times, which dispersed the African population into the forest (Newton , pp. ). Historical documents indicate that the puritan company hierarchy was divided by status di erences in at least three tiers, and each tier had di erential access to agricultural land. Around , governor Robert Hunt received acres of land to farm, whereas men with only moderately high status obtained acres, and men with the lowest rank received acres (Newton , pp. -) (Table ). .
The colonial history of Old Providence was always part of the conflict of colonial powers in the Caribbean Sea. The colonization was conducted by English puritans, and the Spanish crown reacted to this endeavor at least twice between and CE. During these two episodes, Spaniards brought English settlements and the African population on the island under their control (Parsons ; Rowland ). Historic records indicate that a er CE Old Providence was uninhabited for almost sixty years. In the year of CE, there are references of the presence of new Englishmen on the near island of San Andres, but none about permanent occupation on Old Providence (Parsons , p. ). It is only in the year CE that a Spanish o icer documented occupation on the island again. This was represented by only two families from San Andres Island that were living on Old Providence. No nationality or rank of these families is reported, but based on family names of population on San Andres Island, families in Old Providence probably stemmed from the initial Puritan settlements (Ramery de , sheet ). A er that, in CE, Oneill reported about people between free people and slaves living on the island (Oneill , sheet ). Oberto ( , p. ) indicated that in CE, about people lived on Old Providence, of which about two thirds were slaves.

Social Position Description Agricultural Land Received in Hectares
Highest rank English governor . High rank Englishmen . Low rank Englishmen . Lowest rank Runaways ). Areas were classified according to their suitability for the type of communities that settled on the island in the past. The variables used to create the habitat suitability classification are presented in Table . It was assumed that early settlers preferred to locate their dwellings in areas with good agricultural land, abundant terrestrial resources, easy sea access to harbor, and fresh water availability. All variables were weighted as equally important. Each variable was transformed in a four level rank with higher values representing more suitable areas and smaller values representing less suitable areas. Original values of each variable and their assigned values are presented in Table . Drainage basin systems ( Figure ) were used to describe patterns observed in the graphic simulation outputs and input data. .
Marine and terrestrial resources were included in the suitability classification. It was assumed that the most relevant land resource on the island was cultivable land. Areas located in alluvial fans with slight slopes and an average of -rainfall days per year were classified as best for agricultural practices. Other geomorphological units with conditions such as low stability or without watersheds were considered as areas with low suitability to establish dwelling structures or agricultural plots. Black crab (Gecarcinus ruricola) was included as the most important wild terrestrial resource. Today, presence of black crab is concentrated in the western forested areas. Based on a regional study of black crab distribution conducted in (CORALINA-INVEMAR , p. ), the highest value in the suitability rank was assigned to those basins with more than , individuals.  ; Rueda et al. ) suggest that the distribution of fishing practices is homogeneous around the island. Given the small size of the island (17km 2 ), it can be argued that access to marine resources would be determined by access to the sea from the coast. The surface of the island was divided according to the proximity to beach and storm deposits using bu ers of , , and meters wide. These deposits were associated with easy sailing points for small boats, such as those required today to conduct artisanal fishing, and zones closer to these places were classified with higher suitability given the easier access to marine resources they o er.
. Fresh water is paramount for any human population. Although fresh water can be provided through exchange networks or a supply line, these options are limited by logistics. Any sedentarization process supported on non-local fresh water resources will be so limited. This means that permanent human populations seeking self-su iciency require sources of fresh water to ensure their stability; proximity to this resource is therefore an important variable to decide where to locate a residential structure. For this reason, fresh water availability was included in the suitability classification represented by the drainage density of each basin. Drainage density was calculated dividing the sum of the length of the drainages in a basin by the basin area (Horton ). This measure can be used as indicator of surface permeability and proxy for fresh water availability on surface (Table  ). .
The overall suitability classification added the values of suitability ranks for each variable by means of a raster calculation. This classification was originally structured in a x m raster file of the terrestrial area of Old Providence Island. All suitability values were reclassified in four categories: very low suitability ( ), low suitability ( ), high suitability ( ), and very high suitability ( ).

The Model
. The model description follows the ODD protocol (Grimm et al. , ). The description consists of six elements. These elements provide an overview of the model and the general concepts behind the modelâĂŹs design. An animation of one of the runs can be found in the Appendix.
The purpose of the model was to understand how average gregariousness in a socially di erentiated human population changed settlement location in a heterogeneous environment during a colonization process. The model presented did not aim to reproduce historical events. The model simulated settlement patterns under di erent degrees of gregariousness within an enclosed and heterogeneous environment. It was conceived as a tool to model gregarious behavior as a proxy for interdependence and independence.

Entities, state variables and scales
.
The model contains two kind of agents: settlers and runaways. Settlers represent colonial households, and runaways represent escaped slaves. Both settlers and runaways were characterized by the state variables gregariousness, occupation radius, reproduction rate and position (occupied cell). Additionally settlers have the state variable escape rate, and runaways have the state variable fear of settlers (see Table ). Each agent is gregarious only with agents of the same class. Runaways attempt to stay away from settlers.

.
Of all state variables, only position is dynamic: gregariousness, occupation radius, reproduction rate, escape rate and fear of settlers are static variables i.e. they are assigned at agent birth and their value does not change a erwards. Other than gregariousness, which is sampled from a normal distribution, these properties were homogeneous across the whole population. This was a design choice taken in order to emphasize the analysis of the spatial distribution patterns produced by changes in gregariousness. Changing the sampling distribution for any of these attributes can be easily achieved by re-implementing the corresponding sub-procedure in the source code provided.
. The digital landscape in which agents interact is heterogeneous, fixed, and without mobility constraints. It simulates the land extent of Old Providence Island ( . km 2 ) and an adjacent sea area surrounding the Island ( . km 2 ). The total area represents roughly km 2 . Each habitat cell characterizes an area of approximately m 2 . Cell size matches the resolution of the digital elevation model available for the island. Cell size also roughly represents modern rural land property on the island in which % of rural properties are smaller than ha (Rico GarcÃŋa & DurÃąn GarcÃŋa ). Cells are arranged in a D 67 × 101 layout. Habitat cells can be occupied by one agent at a time.
. Landscape heterogeneity is represented by the suitability of cells: the land classification of Old Providence Island introduced above was used to create a four-level rank for cells, representing their suitability for human habitation; inland cells have suitabilities ranging from to , and sea cells have suitability and cannot be inhabited. Suitability here is used only as a representation of the agents' preference for certain terrain, and provides no advantage to the inhabitants. Furthermore, each cell maintains its suitability without changes derived from harvesting activities or environmental hazards. Once again, these factors were excluded from the model in order to focus on the spatial results of interaction between agents with di erent average gregariousness. Each step of the simulation was intended to represent a period of -years.

Process overview and scheduling
. Figure presents the workflow of one simulation step. The following processes take place in each step: • An agent assesses terrain by comparing its position with that of their peers.
• Based on this assessment, the agent chooses between moving or staying.
• If the agent is a settler, they kill all runaways within their occupation radius.
• The agent may then reproduce with a probability given by their reproduction rate.
• If the agent is a settler, it may spawn a runaway through the slave escape submodel (see Submodels).
. This schedule is completed by each agent successively in random order. Agent variables are updated asynchronously while population attributes (means and deviations) are updated globally at the end of each simulation step. Agent attributes are inherited by their o spring. Descendants assess available cells and then settle on an unoccupied cell immediately a er birth. Only runaways die, and only when killed by settlers.

Design Concepts
. ). More complex forms of social interaction such as cooperation could derive from simple rules (Carballo ), like with whom or why you participate as part of a group. In this sense, gregariousness preferences are given in the range between individualism and collectivism, being collectivism a society-level indicator that allows to map communities or groups in which cooperation can emerge. .
Emergence: The main emergent phenomenon is the formation and shaping of communities. Both settlers and runaways change their location and proximity to other agents, and therefore the number and characteristics of their clusters, in response to their interactions.
. Adaptation: Agents identify the optimum cell according to an individual criterion (see objectives) and make a decision about whether to emigrate from their current location or not (see stochasticity). Agents consider suitability of terrain, location of other agents, and their own preference to live near other agents to choose their location. Upon moving to a new cell, settlers exterminate all runaways inside their occupation radius thus shaping the environment to their needs. Conversely, runaways avoid settlers to increase their chance of survival. .

Objectives:
Each agent aims to settle at an optimum location, according to the following criteria: a) Suitability of the terrain (cells) ; b) Neighbors: Agents assess the presence of agents of the same class negatively from the point of view of resource distribution, and positively or negatively from the point of view of companionship, depending on their own gregariousness. With positive gregariousness and as population density increases, the drive towards closeness opposes the impulse to inhabit the best terrain available. Runaways assess the presence of settlers in their neighborhood negatively.
. Prediction: Agents predict the need to split resources (suitability) with other agents inside their occupation radius. Agent predictions do not use memory of previous decisions. Agents predict their gains based on the assumption that other agents maintain the same location. .
Sensing: Agents sense cell suitability and the presence of other agents. They are aware of the satisfaction of agents of their own class with their current locations.
. Interaction: Two agents may not occupy the same cell at the same time. Settlers kill runaways within their occupation radius. Runaways actively avoid settlers. Agents share their location and satisfaction with their current locations.
. Stochasticity: Agent gregariousness is assigned randomly, following a normal distribution. Mean gregariousness and standard deviation are model parameters. The choice to move is a stochastic process, but takes agent's objectives into account. Agents on a disadvantageous position relative to peers of the same class are more likely to move. The choice to reproduce is stochastic: all agents are equally likely to reproduce, with a probability given by their reproduction rate. The escape of a runaway from a settler is probabilistic as well. Other than gregariousness, all random variables are sampled from uniform distributions.
. Observation: Location and gregariousness of each agent were recorded a er each simulation step, along with the total number of agents of each class. . The low suitability of that terrain combined with their own gregariousness keeps settlers away from this area, which allows runaways to thrive -one of the few scenarios where that happens.

Initialization .
Values for the following parameters must be provided in order to run the simulation: • Reproduction rate • Initial population of settlers • Mean gregariousness and standard deviation for settlers and for runaways • Escape rate • Occupation radius • Fear of settlers.
. These parameters were not adjusted with input data; instead, simulations were run with wide ranges of values and the e ect of variations on the parameters was analyzed (See Section ).
. Cells are arranged in a grid and their suitability values initialized according to input from a png file in which suitabilities are encoded as colors. The initial population of settlers spawns in the center of the grid and then their attributes are initialized, one agent at a time, as follows: • Reproduction rate, escape rate, and occupation ratios are assigned the values of the corresponding parameters.
• Gregariousness is assigned randomly, following a normal distribution with the given mean and standard deviation.
• The settler is positioned at its preferred location.
An example of an initialized simulation is shown in Figure (a).

Move .
This submodel is summarized in Figure (a). Agent a assesses the suitability of a cell p as a potential habitat by calculating the sum of suitabilities of all cells within its occupation radius from p, and then dividing the result by the number of agents inside that region. Other agents present will be valued negatively due to the need of sharing produce with them. For this purpose, settlers only take other settlers into account. .
Closeness to other agents may be valued positively or negatively depending on the gregariousness of agents. Only agents of the same class are considered when assessing agent closeness. Last, runaways weight in the presence of settlers within the area. Taking all these factors into account, the overall valuation of the cell is given by: where r(p|a) is the set of cells within the occupation-radius of a from p, A is the number of agents in r(p|a) with which a would share the product of the cells, C is the number of agents of the same as a in r(p|a), S is the number of settlers in r(p|a), and δ a equals if a is a runaway and otherwise.
. Agent a is satisfied with its current position if textitposition − value(a) is greater or equal to the mean position value of agents of its class. In this case, agent a stays in its current position. If textitposition − value(a) is below standard deviation of the mean, agent a is unsatisfied and will move to a new cell. Otherwise, agent a may move to a new cell with probability: where m and σ stand for the mean and standard deviation respectively, calculated on the population of agents of a's class. If the agent moves, it will settle on a random cell with maximum position-value(a).

Slave escape .
A runaway escapes from a settler with a probability equal to the settler's escape rate. The new runaway has a gregariousness sampled from a normal distribution and moves immediately to an optimum location (see Figure  (b)).
(a) (b) Figure : Move (a) and Slave-escape (b) Submodels. Rectangles represent actions, circles represent control switching to another agent.

Simulation Experiments
. Simulation experiments were conducted systematically changing the parameter values within the allowed ranges (see Table ). Only one parameter changed at a time; all the others remained at their default values. A total of replicas were run for each set of values. Each replica consisted of simulation steps.
. The e ect of each parameter on the spatial distribution of settlers and runaways was measured through the Average Nearest Neighbor Ratio (ANNR) (Ebdon ): where d stands for the distance from an agent to its nearest neighbor,d is the average over the corresponding agent set:d and d stands for the expected value of d for a randomly distributed population, which is given by: where A is the area of the minimum enclosing rectangle around the island ( . ). The island's contour was not used to calculate the A because it is very irregular and this a ects ANNR values (Pinder et al. ). That's why we instead considered a quadrilateral area that covered the island while minimizing zones that could not contain population (in this case the sea). Settlers and runaways were considered as two di erent populations, and their ANNRs were calculated separately. .
The ANNRs allow to evaluate changes in spatial clustering for agents of the same class within a fixed study area. If the ANNR is less than , the pattern exhibits clustering. If the index is close to it suggest a random distribution. If the index is greater than , the trend is toward dispersion (Mitchel ). We compared the changes in the spatial distribution of settlers and runaways, simulating long-term changes. We were interested in the variability of ANNR values over time. For each simulation step, the distribution and confidence intervals of ANNR values were calculated using the mean and standard deviation of all ANNR values for each change in parameter values. Parameter values were used as a grouping variable to plot mean ANNR values and standard deviations (Figures , , and ).  Increments in the reproduction rate diminished settler clustering across the digital landscape (Figure (a)). Changes in reproduction rate values for settlers ≤ 1.6% showed an incremental tendency from clustered pattern to random spatial distribution. In runs with reproduction rates ≥ 2% the curve flattens on ANNR= at some point of the simulation, revealing a pseudo-random pattern. For reproduction rates > 2.4% the spatial distribution shi ed toward dispersion in the last steps. This is due to saturation of high-to-mid-suitability terrain enforcing settlement separation outside of the most crowded areas (see Figure ).

Parameter Low Value Used Default Setting High Value Used Increment
. The patterns shown by runaways were di erent (Figure (a)). The spatial distribution of runaways showed a clustered pattern with ANNR values roughly between . and . The runaway population was kept low by settlers' location change: runaway colonies were small and short-lived. In most cases, this pattern continued until the end of the simulation. However, for reproduction rates > 3.2% there was a significant increment in runaway clustering during the final simulation steps. This is due to the increase in settler dispersion and fostered occupation of new areas at high reproduction rates: runaways answered to this situation by fleeing to the lowest-suitability area, which settlers avoided because of its low productivity and isolation (see Figure ). This e ectively created a safe haven for runaways, where an unusually large and tightly clustered community flourished.

Initial population .
The initial population always clustered in the northwest where the largest fertile bay is located (Figure (a)). Historic accounts mention that this place sheltered New Westminster, the first European nucleated settlement. This area is also where the largest nucleated community is located today on the island (Figure ). Agent clustering was sensitive to and increased with the initial population. First settlers and settlers prone to gregariousness were concentrated in the most suitable areas but then population growth resulted in a steady tendency toward dispersion (Figures (b) and (b)). A er the second half of the simulation steps, settler clustering in the simulation with the smallest initial population (n= ) showed considerable spread of its ANNR values. This suggests that the number of possible spatial configurations of agents across the landscape increased. Simulations with larger initial populations showed settler spatial patterns that could be better explained as stochastic processes. Simulations with initial populations ≥ showed ANNR values around , without conspicuous dispersion. In this scenario, neither increase nor decrease in settler clustering occurred but rather low predictability in settler spatial patterns. Runaway spatial patterns did not di er with changes in settler initial population. The ANNR values for runaways were close to , indicating that their spatial configuration probably followed a stochastic process (Figure (b)). These patterns are consistent with those observed for reproduction rate, and are likely caused in the same manner by demographics.

Fear of settlers
. As expected, settlers maintained similar spatial distributions regardless of how much runaways preferred to avoid them (Figure (c)). The ANNR values of settlers increased monotonously during the simulation. Once more, this is explained by population growth enforcing colonization of new areas. Runaway clustering was not sensitive to changes in their preference to avoid cells with settlers either (Figure (c)). Simulation runs with di erent fear of settlers' values maintained a similar tendency and spread of runaway ANNR values, which were more or less constant across the simulation. These values spread roughly between . and , suggesting that runaways had a slight tendency to clustering. This was a consequence of settlers eliminating runaway settlements in each of the new locations they settled, as explained before.

Escape rate .
The increase of escape rate decreases variability in runaway clustering without visible e ects on settler clustering (Figure (e)). For runaways, all escape rate values showed ANNR values below , which could be interpreted as an indication of constant aggregation in space. The lowest positive value ( %) showed more spread of the runaways' ANNR values than other parameter values, at the beginning and at the end of the simulation. The highest escape rates (≥ 40%) created an initial population e ect on runaways that increased the randomness in the spatial distribution of runaways during the first quarter of the simulation steps, a er which the growing settler population overtook runaways (Figure (e)).

Occupation radius .
Simulations were very sensitive to changes in the occupation radius and a high value in this parameter increased settler clustering more than any other parameter included in the model. Low values (between and ) increased the randomness in the spatial locations of settlers, while values ≥ 5 augmented clustering (Figure (e)). This change in regime highlights the increasing importance of gregarious life compared to terrain quality when more opportunities for interaction are available. The spatial location of runaways rather showed a stochastic pattern with increasing spread in ANNR values with respect to the mean (Figure (e)). As before, this is explained by low runaway population due to extermination by settlers. In fact, at minimum occupation radius, where settlers are e ectively blind to the presence of runaways, the runaways population grew and the ANNR curve displayed a di erent pattern, quite similar to the one observed for settlers.

Settler gregariousness .
As with the previous simulations, a trend towards increasing randomness driven by population growth was observed as time progressed. Non-gregarious populations (mean gregariousness ≤ 0) showed almost the same mean ANNR values, falling within each other's standard deviation intervals. On the other hand, increasingly positive gregariousness values (≥ 0.2) visibly increased settler clustering, as expected. Each increment was significant enough to make the mean ANNR values fall outside the standard deviation interval of the next greater adjacent value. Nevertheless, the e ect produced by the population growth through time continued to drive the ANNR values closer toward random patterns (Figure (a)). Runaway clustering was not very sensitive to changes in settler gregariousness. Runaway locations show similar spread in mean and standard deviation ANNR values regardless of changes in settler gregariousness (Figure (b)). Settler and runaway clustering was not sensitive to changes in runaway gregariousness. For both class agents, each parameter showed very similar clustering patterns without significant di erences in the spread of values or the tendency of ANNR values through time (Figures (f), (f)).

Standard deviation gregariousness
. Figures (c) and (d) show that this parameter did not a ect agent clustering for either population.

Discussion
. The di erent scenarios showed that initial settlers start in very suitable areas with a clustered pattern, and new agents move to lower suitability areas o en following a random or evenly-dispersed pattern. Figure (a) shows that during the initialization one major cluster of settlers emerged. Initial settlers found the best spot and thrived there, showing gregarious behavior in resource-rich locations. Figure (b) and Figure (a) suggest that initial gregarious behavior or a large proportion of the population with strong gregarious preferences does not preclude the existence of dispersed agents. On the contrary, simulations showed that on average agents have a tendency to disperse themselves, driven by population growth. Initial settling and environmental heterogeneity limit high suitability location options for subsequent agents (settlers and runaways), increasing their dispersal across the environment. .
Clusters of settlers and runaways can be taken as a representation of local (nucleated) communities (Drennan & Peterson ) and these are intensified by some changes in simulation parameters. Although they did not necessarily located always close to each other, settlers and runaways showed orientation to in-group behavior, which could be interpreted as a sign of collectivism. Initialization always produced one nucleated settlement in the same place where initial settlers of Old Providence started. This place has the largest number of contiguous, highly suitable cells, and for this reason, all simulations started always from a clustered pattern. However, ANNR analysis reveals that the size of this settlement decreases for and negative gregariousness, showing that gregarious drive can play a role in the emergence of large nucleated communities (see Figure (a)). .
Settler clustering was reduced with time in most of the simulations, driven by the population growth of settlers and the restrictions imposed by the initial population and environmental heterogeneity. Runaway clustering showed values in the same intervals through time in all simulations, regardless of parameter changes. Considerable variability in the distribution of runaway clustering, expressed by means of the standard deviation of ANNR values, was observed. This variability could be associated with the dominated role of runaways in the model. Runaways sustained a constant population but could not form permanent groups because they were heavily dominated by settler activities. This suggests that populations structured by dominant roles hinder possibilities of steady interaction between dominated agents. .
Simulations showed that gregarious behavior increased clustering, but the likelihood of interaction between agents showed stronger e ects on agent locations. Settler clustering increased when: a) there was a strong increment in occupation radius (Figure (e) and (b) the initial settler population was very small (Figure (b)). These parameters are modifiers of the number of possible interactions between agents. Large occupation radii trigger more interactions because settlers can intersect with a larger number of other agent occupation radii and cover more settlers with their own occupation radius. A larger occupation radius also increases the number of cells available for each settler which, other things being equal, reduces the disadvantages/benefits of increasing distance between agents and the importance of environmental heterogeneity in ideal location calculations. The combination of these factors averages on a negative or positive e ect on the value of gregariousness. The negative e ect dominated for occupation radius < , causing an increase in clustering. Then, at occupation radius > , the potentiating e ect on gregariousness starts to dominate, causing clustering to increase.
. Runaway clustering, on the other hand, was most a ected by: a) strong increments in occupation radius; and b) high reproduction rate values (Figure (a)). Occupation radius a ected runaways in the same way as settlers, but now another e ect emerged from the survival pressure characteristic of this population: at very low occupation radius settlers' ability to detect runaways becomes hindered, which allows the population to thrive. Similarly, reproduction rate increments increased runaway clustering because new runaways were more likely to emerge from older runaways, but specially because demographic explosion pushed them towards areas less prone to be reached by settlers (See Figure ). This allowed new runaways to cluster near other runaways without being killed by settlers. At the same time, the population growth of settlers created by larger reproduction rates increased settlement dispersal and reduced areas in which runaways could thrive, which in turn decreased the average nearest neighbor distances between runaways that survived a er each simulation step. A small initial population allowed the occupation of highly suitable areas in a slow and orderly process, delaying the emergence of behaviors that locate settlers and runaways in less suitable areas. Simulations thus suggest that the frequency of agent interaction can be driven by the average degree of gregariousness in a structured population but the most conspicuous e ects in clustering were produced by initial conditions and behavioral adaptations that increased the capacity of agents to access more resources and the likelihood of more frequent interaction with other agents. .
The model also suggested that regardless of the agents' preferences, most of the agent population aggregates first within the cells with best suitability and then occupies other less suitable areas. Historic data about Old Providence Island indicate that the colonization process was similar to the pattern in Figure . Historic records describe a mixed settlement pattern configured by a small nucleated settlement, named New Westminster (today Old Town), located in the same area where the largest nucleated community occurred in the model. Dispersed farmsteads and settlements of escaped slaves emerged later on the island. New Westminster only had houses, a church and the governorâĂŹs house (Newton , pp. -). Input data and simulation results suggest that archaeological contexts produced by initial colonial settlements could be located inside the Bowden area. These observations are congruent with theoretical expectations of an Ideal Free Distribution model (Winterhalder et al. ). Simulations also suggest that regardless of the way a human population is structured within a heterogeneous environment, clustering and spatial location are governed by underlying principles such as the likelihood of interaction and the maximization of individual gains, respectively.

Extensions
. The model can be improved in several dimensions. The model assumed a static environment, whereas ecodynamics of human populations and environments create density dependent e ects that are usually expressed in a decreasing relationship (Winterhalder et al. ). In other words, the density of the population can reduce habitat suitability. Dynamic environments with feedback loops that simulate positive and negative e ects of population density is one way to improve the model. Further extensions should explore the e ects of travel costs by means of including variation in the terrain topology. Other model extensions should include stochastic environmental hazards. For example, on Old Providence Island storms and hurricanes occur at irregular intervals and their ecological consequences on settlement nucleation dynamics could be substantial (e.g floods, lost crops). In order to improve the model, extensions can also evaluate the e ects of environmental homogeneity on agent's behavior. The model should be run in environments with less environmental heterogeneity to test its e ects on agent clustering. The modeled scenarios were similar in demographic terms. The e ects of slavery in the model can be expanded by including slave labor and di erent slave trade scenarios as part of the model. The population growth of settlers was constant and this agent class did not die in the model. The model also lacks variables to measure fitness. Agent reproduction can be improved by including death events for settlers and agent fitness measures. The occupation radius is an important parameter in the model and it can be modeled with more variability for each agent class. Di erent occupation radii between or within the same or more agent classes could help to explore complexity behind several levels of hierarchical formations.

Concluding remarks .
This paper describes an exploratory exercise in archaeological agent-based modeling. A simple socio-ecological model yielded results that can be meaningfully compared with historic records and theoretical expectations. This is an encouraging finding. It might be interesting to explore what kind of agent behavior and emergent phenomena occur in larger areas with less environmental heterogeneity, and larger populations with less complete knowledge about their environment merits. Human colonization events limited by technologies of communication and transportation occurred very o en in the past. The archaeological context derived from these processes hold data that could be used to analyze, validate, and plan future long-term human colonization endeavours. One possible future application of this type of model can be found in considering long-term social variables of small hierarchical formations in the future human colonization processes into heterogeneous environments outside planet Earth. The model also shows that even with highly simplistic assumptions, agentbased modeling can be used to assess the e ects of di erent factors in spatial patterns of a colonization process. Di erences in the frequency of face-to-face interaction, central to the modeling process, can resemble cross-cultural variability of individualism versus collectivism found in innumerable studies (Minkov ). The spatial patterns found across space and time in a human community are a proxy for di erent the degrees of gregariousness that a community may show. They capture the spatial aspect of collectivism. Improving the fit of such agent-based models to analyse high level concepts such as gregariousness is still a challenge, but such a broad-brush approach as this one could be a good starting point.