An Agent-Based Model for Simulating Inter- Settlement Trade in Past Societies

Social and computational archaeology focuses largely on the study of past societies and the evolution of human behaviour. At the same time, agent-basedmodels (ABMs) allow the e icientmodeling of human agency, and the quantitative representation and exploration of specific properties and patterns in archaeological information. In this work we put forward a novel agent-based trading model, for simulating the exchange anddistributionof resources across settlements inpast societies. Themodel is part of a broader ABMpopulated with autonomous, utility-seeking agents corresponding to households; with the ability to employ any spatial interaction model of choice. As such, it allows the study of the settlements’ trading ability and power, given their geo-location and their position within the trading network, and the structural properties of the network itself. As a case study we use the Minoan society during the Bronze Age, in the wider area of “Knossos” on the island of Crete, Greece. We instantiate two well-known spatial interaction sub-models, XTENT and Gravity, and conduct a systematic evaluation of the dynamic trading network that is formed over time. Our simulations assess the sustainability of the artificial Minoan society in terms of population size, number and distribution of agent communities, with respect to the available archaeological data and spatial interactionmodel employed; and, further, evaluate the resulting trading network’s structure (centrality, clustering, etc.) and how it a ects inter-settlement organization, providing in the process insights and support for archaeological hypotheses on the settlement organization in place at the time. Our results show that when the trading network is modeled using Gravity, which focuses on the settlements’ “importance” rather than proximity to each other, settlement numbers’ evolution patterns emerge that are similar to the ones that exist in the archaeological record. It can also be inferred by our simulations that a rather dense trading network, without a strict settlement hierarchy, could have emerged during the LateMinoan period, a er the Theran volcanic eruption, a well documented historic catastrophic event. Moreover, it appears that the trading network’s structure and interaction patterns are reversed a er the Theran eruption, when compared to those in e ect in earlier periods.


Introduction
. Social archaeology focuses on the construction of the social space, the exploration of social processes in past societies; and on coming up with models, hypotheses and interpretations on how individuals or groups experience and influence their own societies, and how these constructed narratives of the past can or do influence modern societies. Likewise, the neighbouring discipline of computational archaeology focuses on the utilization of computational and computer-based methods for the study of behavioural evolution, allowing social scientists to model and carry out research in a formal and comprehensible way. Over the past two decades, the social and computational archaeology has utilized agent-based modeling (ABM) for simulating past societies socio-ecological processes based on archaeological information and hypotheses (Axtell et al. ; Janssen . The rest of this paper is structured as follows: Section provides a review of existing examples of archaeological ABMs in the literature. In Section we provide a concise overview of the ABM, by describing its environmental representation; the agents and their interactions; and the agents social organization paradigms. Section presents the theoretical background of the modeling process that was followed for developing the trading network across settlements, based on both the XTENT and Gravity spatial interaction models. There, we also introduce several concepts from network and graph theory required for the analysis of the resulting trading network. Section then presents our specific case study of early Minoan societies located at the wider central area of Knossos in the island of Crete. In addition, we record the empirical evaluation of the various trading models, in terms of potential settlement centralization and organization emerged during the Minoan period, by first detailing the simulation parameters for the various scenarios considered, and then analysing the obtained results. Finally, Section concludes this work, and discusses simulation results and future research directions.

Background and Related Work
. Social archaelogy seeks to understand the social organization of past societies at many di erent points in time (Renfrew & Bahn ). There is a great number of questions that may arise regarding the nature and internal organization of the society under study. For instance, are the main social units, individuals or groups, participate on a more-or-less equal basis, or do prominent di erences in status or rank within the society (or perhaps even di erent social classes) exist? Answering these questions turns out to be harder when exploring prehistoric communities rather than later ones, as written records in early societies are scant. Surely, there are many methods for acquiring information regarding the internal social organization of an early society. Beyond field survey -which aims to discover mainly a presumed settlements hierarchy -one can make use of settlement pattern information, written records, oral tradition and approaches from ethno-archaeology (Renfrew & Bahn ). Clearly, the variety of methods used and the inherent uncertainty of the domain gives rise to a rich space of hypotheses for any given question regarding the social organization of early societies. .
As such, the past few decades has seen archaeology taking an increasing interest in ABM (Doran et al. ; Dean et al. ; Axtell et al. ; Janssen ; Heckbert ; Cockburn et al. ; Crabtree et al. ). Its emerging popularity is due to the ABM's ability to represent individuals and societies, and to encompass uncertainty inherent in archaeological theories. The unpredictability of interaction patterns within a simulated agent society, along with the strong possibility of emergent behaviour, can assist archaeology researchers to gain new insights into existing theories. .
For example, a relatively recent ABM for understanding how prehistoric societies adapted to the American southwest landscape of their era is presented by Janssen ( ). That ABM could explore to some extent how various assumptions concerning social processes a ect the population aggregation and size, and the dispersion of settlements. In that model, interactions (like the sharing of resources among the agents, or the exchange of resources among their settlements) are largely determined by rules that are built in the system. However, results suggest the temporal and spatial population dynamics are a ected by many assumptions of the model; resource dynamics a ect the long-term population levels, whereas climate variability a ects the short-term aggregation levels of the prehistoric populations in American Southwest. .
Likewise, Cockburn et al. ( ) developed an agent-based simulation of specialization in resource production, and a system of barter allowing exchange of specialized products between household agents, in Pueblo societies existing between -CE in southwestern Colorado with the ultimate goal of examining their e ects in these middle-range Neolithic societies. Much of the dynamism of the simulation is provided by annual and spatially specific estimates of potential maize production on this landscape (Kohler ). Agents are responsible for gathering and exchanging resources with other agents. In the case that agents cannot obtain all the resources they need on their own, they are allowed to exchange with other agents. Moreover, if an agent is not performing well at its present location, it will move to a more suitable location in the study area, by evaluating resource productivity of prospective areas. When an agent has excess storage of a resource, it will let other neighbouring agents know; and if they have higher production costs for these resources, they will predict that the agent will be able to help them meet any shortfall for this resource in the upcoming year. Although the authors are less concerned about the realism of some of the assumptions that are hard-coded in their model, they conclude that adding social influence to their reference system with economic specialization leads to a population that is not necessarily larger in size, but one that is more specialized. .
Another recent example of a simulation model integrating cellular automata and a network model of the ancient Maya social-ecological system, is MayaSim (Heckbert ). The purpose of the model is to better understand the complex dynamics of social-ecological systems, and to test quantitative indicators of resilience as predictors of system sustainability or decline. Agents representing Maya settlements (rather than households), develop and expand within a landscape that changes under climate variation and anthropogenic pressure. Agents are utility-based in the sense that they estimate the utility of the various actions at hand. However, they choose actions whose utility has reached some thresholds that are in fact hard-coded by the modeller. Moreover, their utility function is a ected by agent resource exchange that occur between settlement agents, since they are connected via a network of links that represent trade routes. The model was able to reproduce spatial patterns and timelines somewhat analogous to that of the ancient Maya, although this proof-of-concept stage model requires refinement and further archaeological data for better calibrations. .
By contrast, Chliaoutakis & Chalkiadakis ( ) try to assess the influence of di erent social organization behaviours on population growth with respect to their geographic context and available archaeological evidence. Specifically, the ABM's case study is regarding the Minoan civilization, in the "Malia" area at the island of Crete, Greece, during the Bronze Age period. The ABM incorporates several social organization paradigms, giving emphasis to a self-organization social behaviour, where household agents within a settlement continuously reassess their relations with others. An alternative "evolutionary" self-organization social behaviour was also introduced by the same authors, driven by the interactions of strategic agents operating within a given agent community (Chliaoutakis & Chalkiadakis ). They incorporated evolutionary game-theory into the agents social organization process, by simulating repeated "stage games" played by any pair of agents within an agent community, in order to assess the e ects these interactions have on population viability and strategic behaviours that may emerge in the long-term. When agents adopt an "egalitarian-like" social organization paradigm, the emerging development of many "small-size" settlements seems to be the way for survival over time. When the self-organization social paradigm was adopted for determining household agents relations, a "heterarchical" social structure emerges, rather than a clear "hierarchical" one evident in later periods. Simulation results demonstrated that self-organized agent societies appear to be giving rise to larger settlements during their evolution, maintaining larger population sizes than the "egalitarian" distributive one. This fact, is in line with archaeological evidence for larger settlements (towns and palaces) eventually coming to existence during the Middle/Late Minoan period, where a more varied social structure is now suggested (Driessen & Langohr ).
. As a final note, and to the best of our knowledge, the only archaeology-related ABM that utilizes a spatial interaction model, is that of Graham ( ). The ABM simulates movement of travellers (agents) between settlement locations known through archaeological field survey in specific regions of Central Greece during the Geometric period and Central Italy during Protohistory. The author utilizes an entropy-maximizing model, that is, the Gravity spatial interaction model, in order to ultimately rank the settlements by the number of times they emerged as most "important" in the various metrics of the travellers network. Agents in the ABM are only able to travel to settlements around their neighbourhood and only to the most attractive site out of three potential destinations. Although the factual description of the ABM is missing, since the author argues that the mathematics in the ABM are not the most important consideration, but rather the description of how the agents interact, some indicative results are presented and discussed.
A Utility-Based ABM . Against the background presented above, we present here an ABM that serves the following purposes: (a) it can be employed for the study of practically any society of choice in a specific geographic context, and can easily incorporate and assess hypotheses regarding socio-economic processes proposed by archaeologists; and (b) it showcases how MAS-originating concepts, techniques, and algorithms can be employed in archaeologyrelated ABMs. Unlike most existing ABM approaches in archaeology, which employ a simple reflex agent architecture, the ABM here employs utility-based agents that act autonomously towards utility maximization, and can build and maintain complex social structures. Indeed, using ABMs that were built on knowledge derived from archaeological research, but do not attempt to fit their results to a specific material culture, allows for the emergence of dynamics for di erent types of societies in di erent types of landscapes. This aids the generation of new hypotheses, as well as the transfer and re-use of knowledge beyond a specific case study. .
In this work, we adopt and build on top of the work of Chliaoutakis & Chalkiadakis ( ). Agents correspond to households, each containing up to a maximum number of individuals (people). Each household agent resides in a cell within the environmental grid, with the cell potentially shared by a number of agents. Adjacent cells occupied by agents make up a settlement -and there is at least one occupied cell in a settlement. Households are utility-based autonomous agents who can settle (or occasionally resettle) and cultivate in a specific environmental location. They also possess an explicit representation of the environmental grid, and use this to choose the best available migration location they can move to within a pre-specified radius, to improve their utility.

.
The total number of agents in the system changes over time, as individuals belonging to household agents born or die. Agent procreation ability, determining the annual levels of births of household inhabitants, is based on the amount of energy consumed by a household agent x during the year. This in turn depends on the resources harvested, that is, the agent's utility U x , which is a function of population size and available resources at a given location, as we explain below. When household inhabitants exceed a critical number, new households (agent o springs) are created; and when the agent overall utility levels are not high enough to sustain its people, the agent considers migrating to prevent individuals from dying. .
Specifically, at every time step, agent x picks the best action b in the set of actions Actions x at its disposal: The main preoccupation of the agents is to stay alive by acquiring and consuming resources. If an agent fails to acquire enough energy it will eventually die out, since it will stop procreating. Acquiring energy is the only inbuilt goal of the agents. Currently, the only action available for agents to acquire energy is via harvesting the lands. This can be done (a) either at the agent's current location (via employing the cultivation action); or (b) at some other location to which the agent migrates (migration action). Therefore, since there are only two actions to consider, the (expected) utility U x of the agent x can be simply described as follows (assuming the agent cultivates n environmental cells): where R k is a function that describes the agricultural production quantity or reward of a cell k, given its geomorphological characteristics with respect to its location on the grid (land suitability) and the amount of labour applied on the cell by the agents (soil fertility/quality). Equation thus determines that the utility of agent x depends on the expected agricultural production of the cells it cultivates (its total harvested resource amount), as well as the expected utility U x of a new candidate migration location, which in turn depends on the agricultural production of the new position (immediately a er migration). The number of cells n that a given household agent x needs (and is able) to cultivate at a given position, depends on the number of its population (people), and its agricultural regime in use. There are two well-known agricultural practices implemented in our model: "intensive agriculture", where agents intensively cultivate the neighbouring land area leading to greater production per hectare; and "extensive agriculture", where agents "expand" their cultivating area by using more land, but producing less per hectare when compared to the intensive agricultural practice. Thus, an agent may cultivate the land within a specified range from its settled location, while be also able to store any surplus resources to its storage, for up to user-defined number of years.
. Now, agents in the model are able to employ several social organization paradigms, that is, resource distribution schemes, such as: • independent, where agents acquire and consume resources independently • egalitarian-like, where acquired resources are pooled each year and distributed equally among the agents • self-organized, where agents autonomously re-arrange their relations, and hence the underlying "social network structure" describing these relations, without any external control • hierarchical, where agents distribute resources based on "static" relations, following the same rules as those governing the self-organized behaviour but without changing any relation during their evolution . Both in the "self-organized" and "hierarchical" social organization paradigms, agents relations determine the way resources are ultimately distributed among them. .
The self-organization algorithm has two main stages: the task execution and re-allocation mechanism, by which it is decided which agents will receive (energy) resources from others to cover their needs, based on their relations; and the re-organization (decentralized structural adaptation) one, used for re-evaluating and potentially altering their relations at every time step. The evaluation that is performed during the re-organization stage is based on the relative di erence of the overall utility "surplus" exchanged between a pair of agents, in case the relation had been di erent than the current one. Specifically, any interaction between a pair of household agents within a settlement, takes place based on the following relation types: acquaintance, peer or authority (superior -subordinate) related agents, where these relations give rise to a social structure reflecting the flow of resources during exchanges among the agents. The authority relation depicts "superior status" of an agent x over the subordinate agent y in the context of their social organization, reflecting that higher amounts of resources flow from x to y during exchanges than those flowing in the opposite direction; the peer relation holds between agents who are considered more-or-less equal in status (i.e. flows involve resource transfers of almost equal amounts in both directions); while acquainted agents are aware of each other's presence, but have no interaction. Agents use the information about all their past and current year allocations to re-evaluate their relations with others. Thus, agents may improve their performance as a "group" (vitality of the settlement) by modifying the social structure through changes to their relations (re-organization) continuously over time (Chliaoutakis & Chalkiadakis ).

Modeling Trade Across Settlements
. In this work we extend the intra-settlement exchange of resources (household level) to inter-settlement interactions (settlement level). Certainly, in the absence of written records it is not easy to determine what were the mechanisms of trade, or what was the nature of the exchange relationship. However, several formal techniques are available for the study of trade, such as the development of a distribution map for finds or materials, within a specific geographic area (Renfrew & Bahn ). Considering such distribution maps, pondered by fall-o analysis, the quantity of a traded material usually declines as the distance from the source increases. .
For instance, let us consider a "down-the-line" trading system (Renfrew & Bahn ). If one site, e.g. village, receives its supplies of a raw material down a linear trading network from its neighbour site up the line, it may retain a given proportion of the material for its own use, and trade the remainder to its neighbour site down the line. If each village does the same, an exponential fall-o curve will result, as illustrated in Figure . In some cases, however, there are regularities in the way in which the decrease occurs, and this pattern can inform us about the mechanism by which a material reached its destination. As an example, a di erent distribution system, through major and minor sites, would produce a di erent fall-o pattern, in particular, a multi-modal fall-o curve, since lower-order settlements tend to exchange with higher-order centres, even if the latter lies further from the source than an accessible lower-order settlement ( Figure ). We note at this point that in the rest of this paper we shall use the term "settlement" to refer to any site category, such as village, town, etc. .
A possible solution to conceptualize exchange and distribution of resources (flows) between settlements, relies on using a spatial interaction model (Renfrew & Bahn ; Rodrigue et al. ). The basic assumption regarding spatial interaction models is that flows are a function of the attributes W i of the origin location i, and the attributes W j of the destination location j and the "friction" of distance D i,j between the concerned origin and destination locations. The general formulation of the spatial interaction model is as follows (Rodrigue et al. ): In our work here, I i,j represents a measure of "attractiveness" corresponding to the interaction probability or probability of trade between settlements i and j. D i,j is the distance between the settlement locations. Variables W i or W j are used to express a measure of "importance" for settlement i and j, respectively. Attributes o en used to express such variables are socio-economic in nature, such as population or gross domestic product in modern societies. .
Since we are calculating the settlements' interaction probability at any given time-step t during the simulation, we consider the following attributes: • P j,t , defined as the ratio of the population of settlement j (i.e., the total number of its household's inhabitants) with respect to the total population at time t, and • K j,t , defined as the ratio of the number of time steps that settlement j has existed so far up to t.
. Then, at any given time-step t, we define the weight of importance W j,t of a settlement location j as follows: . For example, assume that at time-step t = S i settlements exist in the ABM environmental area, where i = 1, 2 and the total population is inhabitants; and assume that S 1 has a population of inhabitants and a lifetime of years while S 2 has a population of inhabitants and a lifetime of years up to current (annual) time-step t. Then W i,t is calculated as follows: W 1,1000 = P 1,1000 · K 1,1000 = 2880 8000 · 810 1000 = 0.6 · 0.9 = 0.54 W 2,1000 = P 2,1000 · K 2,1000 = 5120 8000 · 360 1000 = 0.8 · 0.6 = 0.48

.
Thus, settlement S 1 has a higher weight (of importance) than settlement S 2 , even though S 2 has an almost double population size than S 1 , due to the higher lifetime of S 1 during the simulation.
. Now, past societies of the first farmers in di erent parts of the world, may be generally described as independent sedentary and relatively egalitarian communities without any strongly centralized organization (Renfrew & Bahn ). Following the development of farming, in many cases, the farming economy underwent a process of intensification, associated with developing exchange. Given this, we make the following assumption: at any given time-step t, each (household) agent within a settlement i is socially contracted as a community member to give away a portion of its stored surplus ps (e.g. % or %) to be communally pooled as the corresponding settlement trading resources N i,t and be traded away by the settlement later on. We note that the percentage of surplus resources that an agent is able to give away is user-defined in the ABM.

.
For instance, if at time t = , settlement S 53 has a i households (agents), where i = 1, 2, 3, and each agent has st i surplus resources in its storage, e.g., st 1 = , st 2 = , st 3 = , while the user-defined percentage of stored surplus to be given away is ps = %, then the settlement's overall trading resources unit N 53,1400 are calculated as follows: Then, settlement i can ultimately trade and exchange resources E i,j,t with settlement j at time-step t, by distributing its trading resources N i,t based on its interaction probability I i,j,t , as follows: To give some intuition on the calculation of E i,j,t let us provide another example; however, in order to not overload notation, we are dropping the t index, when this is not required. Thus, if we consider a set of potentially interacting settlements S i where i = {1, 2, 3, 4, 5} and I i,j is provided by some spatial interaction model (e.g., the XTENT or Gravity used in this work), so that I 1,2 = 0.2, I 1,3 = 0.6, I 1,4 = 0.8, I 1,5 = 0.4 then settlement S 1 will distribute a portion of its trading resources, e.g., N 1 = (in kg) to settlement S 2 , as follows: E 1,2 = I 1,2 · N 1 5 j=1 I 1,j = 0.2 · 200 0.2 + 0.6 + 0.8 + 0.4 = 20 .
As such, S 1 will give away % of its overall trading resources to settlement S 2 , % to settlement S 3 , % to settlement S 4 and % to settlement S 5 -in the event that trade occurs with the corresponding probabilities. Similarly, when the trading process is over, settlement i will proportionally distribute the "public good" payo among its household agents, based on their number of inhabitants.
. We elaborate on the XTENT and Gravity spatial interaction models immediately below.
The XTENT model . The XTENT model asserts some relationship of settlement size and distance, whereby the larger dominates the smaller if the distance between them is su iciently small, whereas the smaller retains autonomy if that distance is large enough (Renfrew & Level ). Thus, it assumes that a large centre will dominate a small one if they are close together; in political terms the smaller site has no independent or autonomous existence. This approach overcomes the limitation of the Thiessen polygons method, where territories are assigned irrespective of the size of the settlement, and where there are no dominant or subordinate settlements, allowing a simple approximation of the political reality and a hypothetical political map to be constructed (Renfrew & Bahn ).
. In our ABM, the "attractiveness" determining the level of trading interaction of settlement i (origin location) with settlement j (destination location) that relies on the XTENT formula, is proportional to the importance of the destination location and declines linearly with their distance, as follows: where β and m are constants used to adjust the required level of the e ect that the weight of importance W j of settlement j and the distance D i,j have on the overall "attraction" between settlements i and j, respectively. Of course, one has to experiment with specific values for β and m to reflect the required attraction between settlements i and j. Intuitively, this "attraction" I i,j corresponds to a (trading) "interaction probability" between i and j. Moreover, in order to turn I i,j into an actual probability of trade between settlements i and j we choose to scale its value to [ ; ] (min-max normalization). .
Given the I i,j 's, one can provide visualization intuitions about settlement main trading territories by coloring each "cell" in the modeling area with the same color of the settlement to which it is mostly attracted to (see Figure in Appendix A).

The Gravity model .
The Gravity model is the most common formulation of the spatial interaction method (Rodrigue et al. ). It is named as such because it uses a similar formulation as the Newton's law of gravity. The "attractiveness" between the locations of origin i and destination j that relies on the Gravity model is proportional to the importance of the destination location and inversely proportional to their respective distance (Hu ): Likewise, one would of course need to experiment with λ in order to e iciently reflect the required (growing) e ect that distance have to the trading probability between settlements i and j. In our simulations experiments and same as with the XTENT model, I i,j is also scaled to [ ; ] (min-max normalization).
. An example of settlement "trading" territories relying on the Gravity model is provided in Appendix A, Figure , by assuming the same settlement locations as in the previous example (cf. Figure , Appendix A) as destination locations, and origin locations to be any landscape cell in the modeling area.

Discussion on the spatial interaction models used .
In the simulation scenarios described later on, we consider two di erent views on the trading probability between settlements; one favouring the distance between settlements rather than its importance, relying on the XTENT model with β = 1.5 and m = 0.005 (Equation ), and another favouring the importance of settlement locations rather than the distance between them, enabled by the Gravity model with λ = 0.2 (Equation ). We will observe these models' e ect on settlement organization and distribution patterns in our simulation results. .
The aim of assigning the specific values of β and m for the XTENT model, and of λ for the Gravity model, is to adequately model the required trade-o between settlements distance and importance for the specific case study's geographic area described later on (maximum distance of about km). To provide an intuition on the two di erent views on the trading probability between settlements, let us assume that the probability distribution of "importance" for a potential destination settlement is as illustrated in Figure (the blue dashed sinusoidal curve). The corresponding probability distribution of interaction of an origin settlement with the respective destination location is then depicted with the red and yellow curve, considering the XTENT and Gravity Figure : Probability distribution of importance of a potential destination settlement location and the corresponding distribution probability of interaction of an origin settlement location, considering the XTENT model with β = 1.5, m = 0.005 and the Gravity model with λ = 0.2. models, respectively. As shown in Figure , the distance between the origin and destination settlements has a greater role when the XTENT method is employed, while it has a lesser impact when the Gravity model is in use.

Graph theory for trading network analysis .
In our ABM, settlements interact with several other settlements, formulating a di erent trading network at every given time-step during the simulation, based on the enabled trading scheme (XTENT or Gravity model). What we need to explore in such a dynamic trading network of settlements, is whether and to what degree some settlements are more important or central than others, based on their trading interactions; and whether settlements tend to create groups characterised by a relatively high density of trading interactions. Thus, in order to better understand and provide insights on the consequence of patterns of interaction between settlements, we adopt in our work some of the main approaches that network and graph theory has developed. We describe these below. .
To begin with, a trading network can naturally be represented by a graph. A graph consists of a set of points and a set of edges or ties connecting pairs of points. In our case, each settlement location in the trading network corresponds to a point in the graph and each trading interaction (link) corresponds to an edge that connects a pair of settlement locations.
. A fundamental measurement concept for the analysis of network graphs is centrality, that can highlight important information about the network organization and its structure (Freeman ). Centrality index describe point locations in terms of how close they are to the "centre" of the network activity. Thus, settlements who have more interaction ties (edges) to other settlements may be in advantaged positions. Because they have many interaction ties, they may have access to more of the exchanged resources over the network as a whole, and hence are less dependent on other settlements (Hanneman & Riddle ).
. Whenever two settlements trade, they are directly connected by an edge, and thus, they are adjacent. The number of other settlements to which a given settlement is adjacent is called the degree of that settlement. A simple and e ective measure of a settlement's centrality is its degree. Since resources can be exchanged in a single edge direction towards another settlement, the temporal trading network of the ABM is represented as a "directed" graph and it is important to distinguish centrality based on in-degree, from centrality based on outdegree. If settlements receive many interaction ties, they can be described as prominent, or having high prestige, since many other settlements seek to direct resources to them, and this may indicate their importance (Hanneman & Riddle ). Settlements with high out-degree centrality are able to distribute resources to many other settlements, or make many other settlements aware of their resource exchange potential, thus being more influential than settlements with low out-degree centrality; although it might matter to which settlement they are distributing resources, this measure does not take that into account (Hanneman & Riddle ). .
Let us now assume that a potential trading network is formulated with n number of settlements S j (network nodes), at a specific time-step during the simulation in our ABM. This snapshot of the trading network can be represented as a directed graph, where numerous trading interactions occur between settlements. The indegree or out-degree centrality index C D (S j ) is the number of incoming or outgoing trading edges, respectively, for a settlement S j (Freeman ): where, tr(S i , S j ) = 1 if and only if S i and S j interact (trade resources) and thus, connected by a tie or edge; and tr(S i , S j ) = 0, otherwise. The magnitude of C D (S j ) for a settlement j partly depends of the size of the trading network on which it is calculated. However, since our trading network is dynamic and constantly changes during its evolution, it is desirable to have a measure that is independent of network size. Thus, we calculate the relative degree centrality C D (S j ) for a settlement j, which is defined as: The e ect of network size has been removed by normalizing with max C D (S j ) = n − 1, since any given settlement S j can at most be adjacent to n − 1 other settlements in the trading network graph. Overall, the degree of a settlement point can be viewed as an index of its potential trading activity. .
Another view of settlement point centrality, within a "directed" network graph, is based on the frequency with which a settlement S k falls between pairs of other settlements on the shortest or "geodesic" paths connecting them, defined as the relative betweenness centrality index (White & Borgatti ): where g i,j is the number of geodesics linking S i and S j , g i,j (S k ) is the number of geodesics linking S i and S j that contain S k and, b i,j is the probability that point S k falls on a randomly selected geodesic linking S i with S j . Similarly to the relative degree centrality C D (S k ) of a settlement S k , the measure is also independent of the dynamic trading network size, since it is normalized by the maximum betweenness centrality of a settlement S k , that is where n O is the number of settlements with outgoing edges, n I the number of settlements with incoming trading links and n S the number of settlements with reciprocated edges (White & Borgatti ). A settlement point in such a position of high relative betweenness centrality can influence other nearby settlements by holding resources in exchange, exhibiting a potential for control of their distribution. It is this potential for control that defines the centrality of these settlements. .
Now, when centrality is applied to the whole trading network graph, such a measure should index the degree to which the centrality of the most central settlement exceeds the centrality of all other settlements, and it is expressed as a ratio of that excess to its maximum possible value for the network graph containing the observed number of settlement points (Freeman ). Thus, the relative degree graph centrality index varies between and , and is defined as follows: where n is the number of settlement points, C D (S i ) is the relative degree centrality defined above for settlement S i , and C D (S * ) is the largest value of C D (S i ) for any settlement in the trading network graph. The maximum possible sum of di erences in settlement relative degree centrality, max , is reduced to n 2 −3n+2 n−1 = n − 2 for the relative degree graph centrality index (cf. (Freeman )).
. Similarly, the relative betweenness graph centrality index varies between and , and is defined as follows: where n is the number of settlement points, C B (S i ) is the relative betweenness centrality for settlement S i and C B (S * ) is the largest value of C B (S i ) for any settlement in the trading network graph. The maximum possible sum of di erences in settlement relative betweenness centrality, that is, max is reduced to n − 1 for the relative degree graph centrality index (White & Borgatti ). .
Then, high relative in-degree or out-degree graph centrality means that there are few settlements of high importance, or highly influential settlements respectively, in the trading network (and thus the most prominent or influential settlement in the network really "stands out", making the value of the numerator in Equation go up). On the other hand, low relative in-degree or out-degree graph centrality means that there are many settlements with a similar level of influence or importance. Accordingly, high relative betweenness graph centrality means that there are few settlements with high potential for control in the trading network, while low relative betweenness graph centrality means that there are many settlements that exhibit a similar potential for control in the trading network. .
To provide visualization intuitions on the relative network graph centrality (Freeman ) we present a snapshot of the trading network developed during a random simulation run. In the example of Figure , each settlement node in the trading network is depicted with a circle, where its size and color represents its relative centrality value [ ; ], with white color corresponding to the minimum value ( ) and black color corresponding to the maximum centrality value ( ). Figure a illustrates a trading network of settlements with high relative graph centrality, while the one in Figure b shows the same network but with low relative graph centrality.
(a) (b) Figure : (a) High and (b) low relative graph centrality indices of a trading network of settlement nodes, represented as circles and trading connections as links between them. Settlement nodes size and color represent their centrality value, from minimum (white) to maximum (black). .
Besides the above relative graph centrality indices that will be used to evaluate the settlement trading network structural evolution, the degree to which settlements in the network graph tend to cluster together is also examined in our work, by calculating the network's average clustering coe icient (Watts & Strogatz ): where n is the number of settlements in the trading network graph and C i is the number of ties between settlement S i 's neighbours, divided by the total number of possible trading edges between its neighbours. C i represents how connected settlement S i neighbours' are. Thus, the network's average clustering coe icientC measures the degree to which settlements tend to cluster together within the trading network.

Case Study: The Minoan Society in Central Crete, Greece
. Several ancient civilizations existed in the Aegean Sea during the Bronze Age, with the island of Crete being associated with the "Minoan" civilization, which came to dominate the islands and the shorelines of the Aegean Sea. A significant shi in the early Minoans human existence and lifestyle was brought when crop farming was first developed. Previous reliance on a nomadic hunter-gatherer way of subsistence, was in time replaced by reliance on the produce of cultivated lands (Hamilakis ). These developments are assumed to have had great impact on the growth of settlements, encouraging the concentration of local population. .
There is not enough information about what kind of relationships existed between the Minoans or how this ancient civilization was organized before the "Post-palatial" (Late Minoan) period. The sophistication of the Minoan culture and its trading capacity is evidenced by the presence of writing (mostly found on various types of administrative clay tablets). The content of the Minoan texts that have been unearthed is predominantly economic (inventories of goods or resources) and religious. Scholars argue that even if relations among (and possibly within) the various towns and cities continued to be contentious and competitive, a common architectural language was beginning to emerge (McEnroe ). This new architectural language marks the beginning of a specifically Minoan identity, which defines a clear indication that each household was not a self-su icient, totally independent production unit, but that it was involved in exchange. Moreover, for the later Neolithic and Early Bronze Age, stylistic and petrographic analyses suggest a low-volume circulation of ceramic vessels, compatible with "gi exchange" economies, over short and long distances between di erent communities within and occasionally beyond the island (Tomkins & Schoep ). This evidence allows us to conceivably model such relations as resource exchanges. .
We note, however, that we do not intend to generally reduce human relations to exchange, as if human ties to society can be imagined in the same terms as a business deal (Graeber ). Nevertheless, even Aristotle was speculating along similar lines in his treatise on Politics. At first, he suggested, families or households must have produced everything they needed for themselves. Gradually, some would presumably have specialized, some growing corn, others making wine, swapping or trading one for the other (Aristotle ; Graeber ). Therefore, although we do not have a clear picture of how human relations (interactions) were actually formed in prehistoric time periods, we need to have a conceivable conceptual model in mind, and that is done with the simplest possible way: to model trading among them as an exchange of resources -thus, giving us the ability to encode the conceptual model as an ABM encompassing various spatial interaction models for the resource exchange process, enabling us to explore a range of the its corresponding trading network structure in turn. .
In addition, archaeologists argue that Minoan palaces are considered to be one of the central factors in bringing about social transformation in the Minoan civilization (Cherry ). In their view, the construction of Minoan palaces came about through a socio-political "quantum leap" from Chiefdom to State. This leap involved also the introduction of writing, the first centrally organized religion (the peak sanctuaries), and the development of social hierarchy and interacting social networks. Moreover, the size of such "grand" public structures at several sites requires both a considerable population and a social cohesion, and it can reasonably be assumed that there were di erent levels of importance, i.e. a hierarchy of sites (Driessen ).
. Furthermore, a series of changes in the Aegean, and in particular in Cretan Bronze Age society, were triggered by the ca. BCE Santorini eruption (Driessen ). These changes would have caused the breakdown of the Minoan system over the course of a few generations, during th c. BCE. Archaeologists hypothesize that the eruption would have initially caused major problems in food production and distribution, undermining central authority and leading to a process of decentralization; this fragmentation would then have led incrementally to internal conflict. Despite the many destructions and abandonments documented, Minoan culture survived. .
Starting from the above archaeological information about the Minoan society, we shall try to assess the resulting trading network structure over time and its e ect on the Minoan society social organization at the community level, providing insights on settlement clustering and organization during the Bronze Age.

Model environment
. The environment is considered to be the geographic area of the wider region of Knossos, located approximately in central part of the island of Crete. As a result, known habitation sites of the Minoan period where identified, categorized and geolocalized, acquired by the "Digital Crete" project . Agents are located within a x km area with one ( ) hectare cell size for the grid space. Moreover, the environment has also associated data layers representing topographical aspects of the model landscape, such as elevation, slope and aquifer locations, contributing indirectly in agent's decision-making process, like where to settle and/or cultivate (Figure ).

Model instantiation .
The estimated per hectare population for an agricultural Minoan settlement during the modeled era ranges from up to (Isaakidou ). In this work, we assume a density coe icient of people per hectare, that is, the maximum number of inhabitants per grid cell (Driessen ). Moreover, the number of household individuals in a given settlement cell is initialized to a random number between and . As a consequence, the maximum number of household agents per settlement's cell is , i.e., divided by the maximum number of individuals per household, that is . Household and settlements number and location are initialized based on archaeological record, i.e., the number of settlements per scenario is set to , which are located at known habitation site locations.

.
Initial cell resources at a given simulation run are based on archaeological estimates on production yield per hectare (ha) pondered by the agricultural regime employed by the agents. As already noted, agents agricultural regimes, can be either "intensive", producing kg/ha or "extensive", leading to a production of kg/ha Figure : Modeling area and its topographical features at central Crete, Greece on an annual basis (Halstead ; Jusseret ). In this work, we assume that household agents employ an intensive agricultural strategy. .
Agent migration radius, that is, the distance that a household agent can migrate to in one time-step is set to the full environmental area (approx. km). An agent may migrate only to a cell where known habitation sites exist, based on the archaeological survey conducted in the specific geographic area. However, we assume a resettlement cost rc for an agent i, which intuitively reflects the decay of potential resources at destination location with increasing distance: where δ is the distance (in km) of the agent to the respective migrating settlement location. The rate parameter of C function above is defined as 0.005 in order to achieve a relatively gradual decay of destination resources for an agent, i.e. model a resettlement cost of about % of agent resources at km away. .
As a final note, we consider a dynamic population growth, based on the amount of resources consumed by a household agent during the year (Chliaoutakis & Chalkiadakis ). This results to a growth rate of about . % when households consume adequate resources, which corresponds to estimated world-wide population growth rates during the Bronze Age (Cowgill ; Jones ).

Simulation scenarios .
We simulate trading across settlements which cultivate the landscape by employing an "intensive farming" agricultural strategy, and deploying the "self-organization" social behaviour, as described in Section . Various scenarios were taken into account for the experimental setup, with di erent parameterization. Specifically, the main simulation scenarios are for our: • two spatial interaction models, the XTENT and Gravity ones, and • two di erent ways to characterize the importance of settlements, one based on Equation , and one based solely on available archaeological data ("site category bias" below) .
We note that a natural disaster sub-model (Chliaoutakis et al. ) is also incorporated in the ABM, in an attempt to provide insights to whether the e ects of the volcanic eruption of Thera (Santorini) a ected the trading network behaviour. The natural disaster sub-model takes e ect at BCE, that is, approximately the date of the eruption estimated by earth scientists (Driessen ). Tsunami and volcanic ash impact on the artificial society and their e ects on agriculture and human life were considered during the conceptualization of the model.

.
Specifically, a meters sea-level rise (including m rise on today's elevation), with inundations of meters inland are assumed for defining tsunami a ected areas on the ABM's environmental grid, rendering associated agricultural fields useless for up to years. We assume that volcanic ash and pumice has the e ect of limiting the production of agricultural fields in the whole model area for up to years. Specifically, ash and pumice are considered to a ect environmental cells inversely linear to elevation, since the volcanic ash layer is smaller at higher elevations-i.e. the impact on agricultural output is % at zero elevation and % at locations with maximum elevation of the area modeled.
. Furthermore, the human impact immediately a er the Theran volcanic eruption, is assumed to consist of % fatalities (a mortality rate due to the event) at the whole environmental area due to one or more earthquakes that the eruption was preceded by, and also due to large amounts of ash and pumice that were emitted. Thus, at the time of the catastrophic event, each individual in our model has a % probability of dying.

.
Simulation results were averaged for each (annual) time-step. Thus, for each scenario, results are averages over simulation runs (repetitions) across a period of , years (cf. Table , Appendix B, for the conventional chronology dates (BCE) of the Minoan period used in our ABM simulation scenarios). Moreover, in all figures below, we depict shaded areas that correspond to % confidence intervals around lines corresponding to agent or network characteristics. In order to assist the reader, in all figures the legends are also ranked in accordance to the relative performance of the corresponding agent or trading network behavioural characteristic. .
The random number generators introduced in parts of the model are obviously "pseudo-random". Thus, via using the same random "seeds", one may introduce the same opportunities for agents in the model simulations, making thereby our simulations reproducible by any interested party. In terms of simulation time, the process can be quite expensive, since a single run (composed of , yearly time steps) takes approximately hours on a single core . GHz computer. However, by employing additional computational power, the simulation process can be sped up significantly; e.g., via allocating a dedicated single-core node of a grid computer to a run, all simulation runs mentioned above can be completed within a few days. Specifically, we utilized the Grid infrastructure of our university, by executing all the above simulation runs on thirty ( ) dual-core ( . GHz) nodes (with GB ram each) in just days (as opposed to months on a single-core computer). .
We now proceed to discuss our findings regarding the trading network analysis performed on our area and era of interest, based on the spatial interaction models enabled and the available archaeological data.

Simulation results
.
We begin with presenting our findings regarding the e ect of the di erent spatial interaction models on household agent population, settlements number, and their size. Simulation results are presented in Figure for both the XTENT and Gravity models, considering a low percentage of stored surplus trading scheme, i.e. ps = %, while agents in the model can settle or migrate only to known archaeological site locations at any specific timestep. The % ps value is in our view a realistic assumption for the age and subsistence regimes studied, given that no sea trade is modeled in this work.
(a) (b) Figure : (a) Number of settlements and (b) settlements size over , yearly time-steps (Minoan period), considering the XTENT and Gravity spatial interaction models. .
When the XTENT spatial interaction model is used, we observe that the number of settlements remains almost constant until the end of the Early Minoan (EM) period, and then gradually increases over time, especially during the Middle Minoan (MM) period and even more a er the volcanic eruption and Late Minoan I (LM I) period (Figure a). The number of agents (households) per settlement also appears to increase until the end of the EM period, and then gradually drops in the MM period. Immediately a er the volcanic eruption, settlement sizes abruptly drop for a few decades, and start again to gradually increase during LM II and LM III periods (Figure b).
. Then, when the Gravity model is employed, we observe a similar behaviour with that of XTENT for settlement numbers and sizes, although the number of settlements is slightly lower than the XTENT model during the EM period, and then slowly increases over time, until the end of the MM period (Figure a). For both spatial interaction models, however, we observe an increase on settlements number and a gradual decline in settlement sizes during the MM period, due to the availability of a lot more known site locations for migration (cf. Figure , Appendix B). We also observe a relatively constant number of settlements a er the volcanic eruption until the end of the LM period, with the XTENT model having higher numbers at about settlements and the Gravity model at about settlements on average. On the other hand, the number of agents per settlement is slowly increasing a er the volcanic eruption until the end of the LM period, with the Gravity model achieving higher numbers of household agents on average than the XTENT model (Figure b).

.
Overall, a higher number of settlements is observed a er the EM period, with an in-parallel decline on the number of households (agents) per settlement. The increasing trend of settlement numbers is in line with the archaeological record, at least until the LM I period, when actual settlement numbers abruptly decline until the beginning of the LM II period (see Figure in Appendix B), and then start to increase again until the LM III period.
(a) (b) Figure : (a) Population and (b) average storage of household agents over , yearly time-steps (Minoan period), considering the XTENT and Gravity spatial interaction models. .
We also report that the overall household agent population is constantly increasing at a dynamic population growth rate from about initial agents to about and agents for the XTENT and Gravity spatial interaction models, respectively, with only an abrupt and short decline immediately a er the volcanic eruption ( Figure a). Furthermore, the stored surplus of the agents is gradually decreasing during the whole simulation period, from about one ton to one half of a ton per household for both the XTENT and Gravity spatial interaction models, with only an abrupt increase immediately a er the volcanic eruption of Thera; and then again gradually decreases until the end of the LM period (Figure b). This "shock" on the average storage of households immediately a er the volcanic eruption appears to ultimately a ect the settlement trading network, since changes in clustering and centralization rates are observed during the LM period, as will be explained later on.

.
Let us now proceed on the study of the structural behaviour of our settlement trading network. In Figure we present the average relative in-degree and out-degree network graph centralities during the , years simulation period. When the XTENT model is employed, the relative in-degree graph centrality gradually drops from about % to % until the end of the EM period (see Figure a) while the relative out-degree graph centrality gradually increases from about % up to % in the same time period (see Figure b). Therea er the relative out-degree graph centrality gradually declines to about % until the end of the MM period, abruptly declines immediately a er the volcanic eruption to about % and then again increases to up to % until the end of the LM period. Relative in-degree graph centrality is kept almost constant to about % until the end of the LM period, with an abrupt and short decline immediately a er the volcanic eruption. . Low relative in-degree graph centrality rates observed during the EM and MM periods (under XTENT) suggest that there are no clearly "prominent" settlements, meaning that, there are no central attractors considering the other settlements in the trading network. On the other hand, the in-parallel high relative out-degree graph centrality rates during the same period, indicate that there are a few settlements that are considered influential in terms of resource distribution. Therefore, one could assume that a settlement organization of distributing resources by these influential settlements in the trading network is implied, at least before the volcanic eruption of Thera or the LM period.
. Using the Gravity model, the relative in-degree graph centrality gradually increases from about % to % until the end of LM period, however, with an abrupt fall and rise immediately a er the volcanic eruption (Figure a). By contrast, the relative out-degree graph centrality slowly decreases from about % to % during the whole period, with an abrupt decline immediately a er the volcanic eruption (Figure b). These high relative in-degree graph centrality rates (under Gravity) suggest that there are only a few "prominent" settlements in the network, implying the possibility of a settlement hierarchy where resources are traded towards these important settlements by other settlements in the trading network. Notice however that this assumes an "attractiveness" of the sites given their W i importance defined via Equation , and not the known category of the archaeological sites.
In the next section, we see that the "conclusions" obtained with the Gravity model are quite di erent when the real sites' category is taken into account; and that in that case they are more in agreement with those of the XTENT model. .
Moreover, the relative graph centrality based on betweenness is considerably low regarding both XTENT and Gravity models, as presented in Figure a. This means that most of the trading connections can be made in the trading network without the aid of an intermediary settlement. Thus, there do not appear to exist settlements with much potential of controlling the inter-settlement trade. As such, there is a need to further study if there are other group formation phenomena at work, which need to be captured. .
Studying the average clustering coe icient of the trading network graph (Figure b), we observe that when the Gravity model is employed, it is relatively low (below %) until the beginning of LM period, while it is relatively high a er the volcanic eruption (more than %) until the end of the LM period. When the XTENT model is employed for the trading process, we observe that the average clustering coe icient of the network graph gradually declines from about % to % until the end of middle EM period; however, it then gradually increases to about % until the end of the LM period, with an abrupt and short fall immediately a er the volcanic eruption.
(a) (b) Figure : (a) Relative betweenness graph centrality and (b) average clustering coe icient of the trading network over , yearly time-steps (Minoan period), considering the XTENT and Gravity spatial interaction models.
. Thus, for both the XTENT and Gravity models, the observed settlement trading clustering behaviour a er the volcanic eruption until the end of LM period, implies a more dense trading activity between settlements at the time, raising the possibility of more settlement clusters in the trading network. Assuming that such settlement clusters were around large towns, cities, or palaces, this trading network clustering behaviour has a correspondence to the archaeological record ( Figure , Appendix B): just two cities are known to have existed during the EM period (Archanes and Knossos), while several large towns, cities and palaces were flourishing in the area during the MM and LM periods (Knossos, Malia, Archanes and Galatas). .
Finally, for interest, we also conducted additional experiments considering the same simulation scenarios, however, with a higher percentage of stored surplus trading scheme, i.e. ps = %. Simulation results exhibit similar behaviour with no remarkable di erences, besides the average storage per household agent, where even lower amounts of resources stored are observed for the scenario of trading a higher portion of stored surplus. Corresponding results figures are presented in Appendix C, since their behaviour is entirely similar with simulation scenarios considering a lower percentage of stored surplus trading scheme. This similarity in the trading behaviour observed in the results where ps = % is justified, since the trading network structure naturally takes into account only the number and density of trading interactions between settlements, and not the volume of resources exchanged within the trading network as well. .
In all of the above simulation scenarios, we used attributes relating to settlement's population and lifetime during the simulation period for calculating the weight of the importance W i of a settlement i, given Equation .
In the following simulation scenarios, we fix the W i values with known archaeological site categories. This will enable us compare the settlements trading network and the organization structure developed, based on archaeological estimates on settlement types, with the one autonomously developed during the simulations described above.

Site category bias .
Let us first assume a simple, broad classification of settlement types rather than specific site categories, which corresponds roughly to the site hierarchy put forward by Driessen ( ), based on archaeological estimates: village (or settlement or hamlet), corresponding to less than . ha catchment size, hosting fewer than households / inhabitants on average; city (or large town or town), corresponding to less than ha catchment size, having fewer than households / inhabitants on average; and palace (or capital town), corresponding to greater than ha catchment size, with more than households / inhabitants on average. Based on this classification of settlement types, instead of using Equation , we express W i of any settlement point location i as a weight in [ ; ], by mapping the corresponding known archaeological site type as follows: • W i = 0.5 when the corresponding archaeological site category is a village, • W i = 0.7 when the corresponding archaeological site category is a city, and • W i = 0.9 when the corresponding archaeological site category is a palace .
As such, the "attractiveness" or the probability of trade for any settlement in the trading network, is biased by the corresponding known archaeological site category. Thus, in the following simulation scenarios, settlement importance is based on archaeological evidence on the settlement type at any given time-step and geographic location. The rest of the experimental setup is exactly the same as the simulation scenarios discussed in the previous section. .
To begin with, simulation results on agent settlements number and size are presented in Figure for both the XTENT and Gravity models. We observe that the number of settlements remains relatively constant until the end of the EM period, similarly to the previous scenarios, where settlement importance was calculated by its own dynamic characteristics, i.e., population.
(a) (b) Figure : (a) Number of settlements and (b) settlements size over , yearly time-steps (Minoan period), considering known archaeological site categories for both the XTENT and Gravity spatial interaction models. .
Regarding the XTENT model, we observe a similar behaviour with scenarios not being biased by site categories, where a gradual increase of settlement numbers over time is noticed, especially during the MM period and even more a er the volcanic eruption and LM I period (Figure a). Similarly, the number of agents per settlement increases until the end of the EM period, and then declines during the MM period. This is due to the high migration rates (because of population growth) observed to more (known) settlement locations available during that period. Moreover, settlement sizes abruptly drop immediately a er the volcanic eruption, however, then gradually increase until the end of LM period (Figure b).
. When the Gravity model is employed, we observe a similar behaviour with the XTENT model in settlement numbers and sizes, although the number of settlements slightly declines at the end of the MM period, and drops further immediately a er the volcanic eruption; and then remains relatively constant until the end of the LM period ( Figure a). Thus, in contrast to the previous scenarios, where no bias by known archaeological site categories was introduced, an entirely di erent behaviour is now observed. That is, a significant di erence in settlement numbers is observed, growing up to about settlements during the end of the MM period, and holding up to about settlements until the end of the LM period, while just a number of about and was observed in the previous scenarios (cf. Figure a). We note that this trend in settlement numbers is surprisingly very similar to the one that exists in the archaeological record for the specific environmental area during the whole Minoan period (cf. Figure , Appendix B), with the only di erence being a substantial decline reported at the end of LM I period in the archaeological record -and which was due to unknown "external" events. Higher values in settlement numbers exist in the archaeological record, suggesting that a higher population growth rate (> . %) probably should have been used during our simulations (we chose to follow Cowgill ( )). .
On the other hand, the numbers of agents per settlement tends to increase until the end of the EM period, and then abruptly declines at the beginning of the MM period from about to households and further decline during the MM III period down to . The number of households per settlement, however, is slowly increasing a er the volcanic eruption until the end of the LM period, with the Gravity model not being able to achieve higher numbers of household agents per settlement on average than the XTENT model (Figure b). .
We also report that the overall number of households (i.e., the agent population) is constantly increasing during the whole time period, same as in the scenarios without bias from known archaeological site categories, being able to even achieve higher population sizes, with only an abrupt and short decline immediately a er the volcanic eruption (cf. Figure in Appendix C). We note that we observe higher population sizes a er the Theran eruption when site category bias is used and the trading network is simulated via the Gravity model (as opposed to XTENT).
(a) (b) Figure : (a) Relative in-degree and (b) relative out-degree graph centrality indexes of the trading network over , yearly time-steps (Minoan period), considering known archaeological site categories for both the XTENT and Gravity spatial interaction models. .
Regarding the structural behaviour of the settlement trading network, the relative in-degree and out-degree graph centralities are presented in Figure . The XTENT model exhibits a very similar behaviour to the one without known site types bias (cg. Figure ). Interestingly, the Gravity model is now showing a similar behaviour to the XTENT model, that is, it exhibits lower rates of in-degree and higher rates of relative out-degree centrality. The low relative in-degree graph centrality rates during the EM and MM periods, imply that there are no "prominent" settlements. By contrast, the high relative in-degree graph centrality rates observed in the trading network a er the volcanic eruption and during the LM period, suggest that there are certain "prominent" settlements in the trading network. On the other hand, the low relative out-degree graph centrality rates during the LM period, indicate that there are many settlements with a similar degree of "influence" in terms of resource distribution. Therefore, one could assume that a settlement hierarchy where resources are traded towards the (few) most important settlements in the trading network is implied during the LM period. .
Moreover, the relative betweenness network centrality is low for both XTENT and Gravity models, as presented in Figure a, even lower than scenarios without site category bias (cf. Figure a), suggesting even less potential for control on the flow of resources traded between settlements. However, there is a structural basis for assuming that certain settlements with the highest relative betweenness centrality in the society are "di erent" from the other settlements in the area, at least during the EM and MM period. .
Regarding the average clustering coe icient of the trading network graph (Figure b), we observe that the Gravity model has again a similar behaviour to the XTENT model, that is, it gradually declines from about % to (a) (b) Figure : (a) Relative betweenness graph centrality and (b) average clustering coe icient of the trading network over , yearly time-steps (Minoan period), considering known archaeological site categories for both the XTENT and Gravity spatial interaction models.
% until the end of middle EM period, and gradually increases to more than % until the end of the LM period, with an abrupt and short fall immediately a er the volcanic eruption. The low clusterization thus observed in the trading network until the end of the EM period may suggest that the trading network connections are losing density until the end of the EM period. The network's clusterization appears to be recovered in the MM period, and even more in the LM period, indicating the possibility of more dense settlement clusters in the trading network, where resources are traded towards the few most important settlements within these clusters (those with high relative in-degree graph centrality). There seems to be a correspondence with the archaeological record, enhancing such a possibility -since several towns, cities or palaces are recorded during the MM and LM period, while just a two towns exist during the EM period, as previously noted. .
Concluding this section, we remind here the reader that the Gravity model is able to better capture the trend in settlement numbers that exist in archaeological record. This is a reason to believe that, in this case, Gravity allows us to better interpret the structure and dynamics of the formed trading network. The "unchanged" behaviour of the XTENT model is justified, since it favours the distance between settlements rather than their importance. Thus, it should be used in cases where settlements importance is not known, or cannot be properly modeled.

Conclusions and Future Work
. In this work, we presented an archaeology-related ABM to simulate trading interactions across settlements of an artificial past society. The ABM was formally built upon MAS-originating concepts, adopting a utility-based agent design. We model inter-community trading interactions by incorporating a trading sub-model, employing two well-known spatial interaction models, XTENT and Gravity. The simulations' aim was to assess the sustainability of the artificial society in terms of population size, number and distribution of agent communities with respect to both spatial interaction models; to analyse the resulting trading network structure during its evolution over time; and to draw interesting conclusions (or, rather, sketch out interesting hypotheses) about the settlements' hierarchy, via annotating our results with the archaeological record. .
As a case study we considered the Bronze Age Minoan civilization and as the ABM's environmental area we considered the geographic area of the wider region around Knossos, located in the central part of the island of Crete, Greece. Simulation results show that when settlements' importance is known or properly inferred (based on archaeological data or evidence), modeling a trading network relying on the Gravity model can produce settlement patterns similar to the one that exist in archaeological record for the area under study (see Figure a).
On the other hand, if solely settlement locations are known, then the XTENT model can produce acceptable results on simulating the trading activity between them.
. Simulation results show that, when the importance of settlements is weighted based on the type of known archaeological habitation sites, high relative out-degree centrality rates are observed in the trading network, along with low clustering during the end of the EM period. This observation suggests that a small number of influential centres might have existed, linked to a settlement hierarchy where resources are distributed by these centres to other settlements in the network -but with no clearly prominent centres to which resources are directed. Interestingly, a er the catastrophic event of the volcanic eruption of Thera and during the LM period, the trading network connections are becoming much denser, and resources are now being distributed towards only a few prominent settlements in the network. We note that these results are in line with archaeological theories suggesting a er the MM period the actual settlements hierarchy was transformed, with subsequent radical changes in their trading network, a ected by settlement numbers and sizes as well as natural disaster events (as also indicated by Figures and in our work here). Specifically, archaeologists argue that independent political units and centres of the EM and early MM period, were incorporated into a "Knossian" state during late MM and early LM periods by being demoted to secondary centres while others were promoted from tertiary to secondary centres (Driessen ). Thus, large and comparatively well integrated polities that existed until the end of the MM period in Central Crete were incorporated into a larger political framework and a territorial state headed by Knossos (Driessen ). Given the above simulation results, our ABM appears to be able to provide support for those theories to some extent. .
In terms of future work, we need to run more scenarios with a variety of initialization setups. In addition, we intend to equip the ABM with additional modules (vegetation data, geological information, other archaeological evidence or scenarios of interest) and further types of utility-generating activities (fishing, hunting, animal husbandry, cra ing), to allow the further integration of qualitative and social factors into our model (Roux ). Specifically, we intend to run simulation scenarios where additional resources during the trading process are available from external sources, by building a more elaborate trading model that will also include maritime trade (and its e ects in coastal settlements) or trade based on the proximity to religious centres or peak sanctuaries; of course, these components will be also a ected di erently by the natural disaster sub-model, when it is enabled (i.e. the volcanic eruption of Thera, in this case study). We also intend to represent the dynamic trading network as a "weighted" directed network graph, in order to take into account the amount or volume of resources exchanged during trade, rather than solely the number or density of trading interactions between settlements in the network. Finally, we intend to run simulation scenarios on (artificial) past societies in different geographical space and time, where su icient archaeological data is available for testing and assessing ABM results with respect to related archaeological hypotheses regarding their social organization.
Households were considered to be the main social unit of production during the period and in the society under study (Whitelaw ).
For details on the definition of R you may refer to Chliaoutakis & Chalkiadakis ( ).
More accurately socio-economic social organization paradigms.
The distance factor D i,j is measured as the Euclidean (linear) distance for simplicity. This distance can be alternatively measured as the Least Cost Path between two settlement locations, considering slope and elevation as cost surfaces, however, with significantly higher computational cost.
See http://digitalcrete.ims.forth.gr Short "jumps" observed in the figures immediately a er the volcanic eruption are not real, but a result of the Savitzky-Golay smoothing filter applied on the data (Savitzky & Golay ). The filter increases the precision of the data without distorting their tendency, by fitting successive sub-sets of adjacent data points with a low-degree polynomial with the method of linear least squares.
We remind the reader that all potential settlement locations correspond to actual settlement sites.
Archaeologists assume that a wave of fire destructions a ected Cretan settlements during and at the end of LM IB, that have variously been attributed to internal revolt, Mycenaean invasion, or to a major natural disaster involving earthquakes (Driessen ).

Appendix A: Visualizing settlements attraction areas
In the following two figures, each settlement is depicted with a unique coloured circle, where its size represents its importance with respect to its type (village, town or city in this example), while its main trading territory, i.e., landscape cells that are mostly attracted to it, is depicted with the same color. Figure : Visualization of "territories" of di erent settlements (of type village, town or city) within the modeling area, considering the XTENT spatial interaction model, considering β = 1.5 and m = 0.005. In particular, settlement (of type village), located near the centre of the modeling area, will most probably trade with settlement (city) or even settlement (village), since it is attracted to settlements that are relatively close in range, undervaluing the importance of settlements that are further away. Figure : Visualization of "territories" of di erent settlements (of type village, town, city) within the modeling area, considering the Gravity spatial interaction model, considering λ = 0.2. Now, settlement (village) will most probably trade with settlement (city) or even settlement (town), since it is attracted with settlements of high importance, that is of type city or town, despite its distance from them.   questions: How a Minoan settlement trading network would look like in the wider area of Knossos, and how does the trading network's dynamics evolve in response to environmental changes and the anticipated social crisis caused by the volcanic eruption of Thera island? For this purpose, we consider two di erent approaches for modeling trading interaction probability, corresponding to "attractiveness" among settlements, based on the distance between them and their "importance". As a case study, we consider the "Minoan" civilization and provide in the process insights and potentially support for archaeological hypotheses regarding the settlement organization in the wider area of Knossos in Central Crete, Greece, during the Bronze Age.
Simulation results show that modeling a trading network by favouring settlements importance rather than the distance between them, as with the Gravity model, can produce settlement patterns similar to the one that exist in archaeological record for the area under study. However, such a spatial interaction model should be deployed only when settlements importance is known or properly inferred. In such case, the trading network structure also suggests that a small number of influential centres could have existed until the Early Minoan period. This fact may suggest a settlement hierarchy where resources are distributed by these influential settlements to others in the network, however, with no clearly prominent centres to which resources are directed. In addition, a er the volcanic eruption of Thera and during the Late Minoan period, the trading network connections are becoming much denser, with resources being distributed towards only a few settlements in the network.
Patterns: Archaeological estimates for the era and space under study are provided as input in the model, as well as archaeological survey data; in particular, known habitation site locations and types, such as village, city or palace. In our simulation experiments, agents in our environmental grid are allowed to settle or migrate only to cells that actual archaeological habitation sites existed at the specific time (step). Therefore, the patterns observed on settlement numbers during our simulations against the ones observed in the archaeological record, allow us to evaluate the suitability of two di erent spatial-interaction models, that we used for modeling the interaction probability (attractiveness) among settlements.

Entities, state variables and scales
Entities: The are three main entities that are included in the model: households, settlements and environmental cells. Agents correspond to households, each containing up to a maximum number of individuals (household population). Each household agent resides in a cell within the environmental grid, with the cell potentially shared by a number of agents. Adjacent cells occupied by agents make up a settlement -and there is at least one occupied cell in a settlement. Each agent cultivates a number of cells located next to the settlement. The number of cells that a given household needs (and is able) to cultivate at a given location, depends on the number of its individuals, and its agricultural regime, i.e., "intensive" or "extensive" agriculture.

State variables:
The main state variables for a household, besides the current number of its individuals and the settlement that the household belongs to, are: current utility, minimum resources needed, surplus of resources stored, a list of its current cultivating patches, a boolean variable (true/false) for considering migrating to another location, the number of cells that needs to cultivate, the number of years (annual time-steps) that the agent exists within the settlement and in the environment, and others.
The main state variables for a settlement are: a list of the household agents within the settlement; the total number of households and its population (number of people); the weight of its importance; the total amount of resources received by other settlements during trade; the number of years that the settlement exists in the environment; the in-degree, out-degree and betweenness centrality values; and the clustering coe icient regarding the settlement's position on the trading network.
State variables for environmental cells include static variables, such as elevation and slope data derived from a today's Digital Elevation Model (DEM) of the modeling area; the land-productivity for the patch (amount of resources) that is based on settlement's population; the land-suitability that is based on the slope value of the corresponding cell; the utility of an aquifer for the respective cell; the agricultural production (agent reward) of the cell that is based on land suitability and productivity; the initial resources of the cell; and several other parameters denoting if the cell is also an aquifer location, or if a known archaeological habitation site existed in the current time-step (along with site type and derived importance), or if it is currently cultivated, and others.
Scales: Household agents and resources are located within a x km area, with a x m ( ha) cell size for the grid space. Thus, the landscape consists K cells. The time slot investigated is , years (ca. -BCE), with annual time-steps.

Process overview and scheduling
Processes: At every time-step, corresponding to a period of one year, household agents first harvest resources located in nearby cells (corresponding to the fields they are cultivating). They then check whether their harvest satisfies their minimum perceived needs. In such a case, potential surplus is added to any stored resource quantities for a specific number of years (user defined parameter). If not, they might ask others for help (depending on the social organization behaviour in e ect), or they might even eventually consider moving to another location or settlement, depending on the expected utility for these actions. Moreover, we assume that each (household) agent within a settlement is socially contracted as a community member to give away a portion of its stored surplus (e.g., %) to be communally pooled as the corresponding settlement trading resources and be traded away by the settlement later on. Then, the settlement can trade and exchange resources with any other settlement based on the probability of interaction or "attractiveness" with the other settlement, thus ultimately transforming the trading network structure. Furthermore, population size a ects the land productivity, in two ways: positively, since the continuous occupation or cultivation of an area by a large populace leads to experience and subsequent higher crop yield; and negatively, since it also leads to overexploitation of resources, thus,induce a lower crop yield. Population levels at a given area are a ected by migration, as well as natural population change by birth and death of household population. Lower amount of resources reduces birth rate and thus leads to a reduced population size and threatens household agents with extinction. We also consider a natural disaster that takes e ect at a specific year ( BCE, estimated by earth scientists) that is, approximately the date of the volcanic eruption of the Thera (the nearby island of Santorini) in an attempt to provide insights to whether the natural disaster a ected the trading network behaviour and how.

Schedule:
The model is developed to cover the needs of the potential organization of an artificial Minoan society, both in the household and settlement level, based on archaeological evidence and data. At every timestep, the following processes take place: first, we get and set known archaeological habitation sites locations and types information on environmental cells from the associated archaeological sites dataset; the respective information need to be available to agents that consider migrating (process) to another location. Then, environmental updates regarding resources and potential cultivation area for the households agents are scheduled. Cultivation (and harvesting) process follows, where the agricultural production, that is agents reward, for each cultivated field is calculated, based on land suitability and productivity. Household agents expected utility is calculated from all their actions (currently cultivation and migration), immediately a er considering the social organization paradigm deployed within the settlement. A erwards, the trading process (resources exchange) across settlements in the model environment is scheduled, based on the spatial interaction model employed (XTENT or Gravity), thus further updating household agents utility along with the trading network structure. Then, updates regarding agents storage, utility, and their associated characteristics are scheduled, followed by updates on population size, a ected by a dynamic population growth rate of . %, when households consume adequate resources. Moreover, population size is further updated along with available agricultural resources, if the natural disaster sub-model is enabled, by simulating the e ects on agriculture and human life of the volcanic eruption of Thera, based on archaeological evidence and estimates. Finally, observer actions, that is plotting graphs and recording necessary variables to output files, take place at the end of every time step.

Design concepts
Our model employs a utility-based agent design that act autonomously towards utility maximization, that can also build and maintain complex social structures. The utility-based agent architecture, introduced by Multi-Agent Systems, can provide a more general performance measure of the agent "well-being", since it allows a comparison of the di erent agent actions-states to exactly how well the agent is in any given time-step (rational utility maximizer). The model also includes several other design concepts:

Emergence:
The main outcomes of the model are trading interaction patterns. In particular, how dense the whole settlement trading network is and if it exhibits an accumulative or distributive behaviour in terms of resources towards a few important or influential settlements, respectively. These outcomes emerge from archaeological data introduced in the model, i.e., the spatial distribution of available known archaeological habitation site locations where agents are only allowed to settle or migrate to and, to a lesser extent, from the way settlements importance is defined and modeled, a ecting settlements interaction probability.
Adaptation: Household agents exhibit an adaptive behaviour: deciding whether or not to migrate to another location and, if so, selecting a new location based on the agent's expected utility for the new location; and also a structural adaptation behaviour when the self-organization social paradigm is employed.
The decision of whether to move depends on satisfying the only in-built goal of the household agent, that is, to stay alive; If an agent does not receive the minimum level of resources it requires (utility threshold) for a specified number of years in a row (and the storage is empty), then it considers migrating to another location (or settlement).
When the self-organization social organization paradigm is employed by the households within a settlement, they interact with one another ("gi " exchange) for the proper allocation of resources. This is achieved by adapting the social structure through changes to their relations (re-organization) continuously over time, thus improving their performance as a "group" (vitality of the settlement).
Sensing: Household agents are assumed to know their own attributes, such as their current utility, surplus stored, etc., and its current relation with other households in the same settlement. Agents have also an explicit knowledge of the environment within a radius that is user-defined (model parameter), although we assume a full environmental view for the agents due to the small area modeled (maximum km).
Interaction: There are two main types of interactions in the model, agent -environment and agent -agent interactions. The former describes the interaction of household agents with their environment, while the later describes their interaction with other households (or settlements) in the modeling area.
The interactions related to agent -environment actions that are currently implemented in our model, are cultivation and migration actions. Both actions depend on the agricultural production of an environmental cell that is a ected by land suitability, expressed as a decay with increasing slope, and land productivity representing soil fertility and a ected by current population size of the settlement (and inspired by the logistic map equation). We note here, that household agents are making a decisions for cultivating a specific number of cells, based on their number of individuals and the agricultural practice in e ect, thus allowing them to achieve an acceptable level of resources instead of maximizing it. With this assumption we are able achieve an appropriate trade-o against utility maximization consideration introduced by the agent design. Likewise, a resettlement cost is also introduced during migration decision-making, that is defined as a decay of expected resources at migrate location with increasing distance. Re-settlement cost is intuitively reflecting costs associated with the relocation of the household agent.
The agent -agent interactions are structured in two levels, intra-settlement interactions that are defined based on the social organization paradigm employed by the household agents within a settlement, and inter-settlement interactions, defined by the spatial interaction model deployed.
Stochasticity: Stochasticity is used in two ways. First, the model can be initialized stochastically, regarding the initial number of settlements in the model, initial number of households, number of household individuals, and amount of resources in each cell. These initialization methods are stochastic so that each simulation run for each scenario produces di erent results. In this work, in particular, only the number of settlements is not stochastically initialized, but rather based on actual number and locations of known archaeological habitation sites at the beginning of our simulation (ca. BCE). Moreover, there is also an inherent stochasticity introduced by the (NetLogo) modeling environment regarding the way that agents are executing related processes in the model; households agents or settlements are executing each and every process in a random order.
Second, population dynamics are stochastically defined and based on the dynamic population growth rate (model parameter) that is a ected by the amount of resources (thus utility) consumed by the household agent during each year. In addition, settlements are trading resources with one another based on their "attractiveness" (interaction probability) that is a ected by the distance between them and their (weight of) importance, and thus the spatial interaction model in e ect. Collectives: This model includes two kinds of collectives: households that are representing groups of inhabitants, and settlements representing groups of households. The collectives are represented as specific entities with their own state variables and behaviours, as defined above at Entities, State Variables, and Scales. These collectives are included in the model because are essential in modeling an artificial human society from the bottom up.
Observation: The model was tested and analysed by recording all entities and model variables, and observe and validate their values process by process. Moreover, at the graphical user interface (GUI) of our model some of these variables are traced at every time-step for quick evaluation, such as the number of households and settlements, population size, average utility and storage (and per settlement), agent migrations and average migrations per (annual) time step, relative graph centrality (in-degree, out-degree, betweenness) and centralization (clustering coe icient), and others. In any case, all non-static variables are recorded and provided in output files for further statistical analysis.

Initialization
The number of households and population size for each household agent are randomly initialized from a sample of a uniform distribution, depending on the maximum number of inhabitants per cell (ha) and individuals per household agent. The amount of resources for each environmental are also randomly initialized from a sample of a normal distribution. In this work, initialization of state variables is intended to be case-specific, that is, an artificial Minoan society at a specific geographical territory of the island of Crete, Greece (the wider area between the archaeological sites of Knossos and Malia). Therefore, the number of settlements is initialized to twenty one ( ) placed at specific environmental cells, that correspond to specific geographical locations of known archaeological habitation sites locations, that are provided as input spatial data file (shapefile).
The maximum number of inhabitants per cell (ha) is set to people and the maximum number of individuals per household is set to . Regarding the environmental resources, the minimum amount of resources required per individual per year was set to kg, while households employ an intensive agricultural practice leading to a production of maximum , kg /ha. All these parameters are based on archaeological estimates for the period and area under study (Halstead ; Isaakidou ; Jusseret ). Moreover, environmental cells state variables, such as elevation and (derived) slope (raster DEM), or aquifer locations (shapefile) are also provided as input spatial data files to the model.
Various scenarios are taken into account for the experimental setup, with di erent parameterization regarding the two spatial interaction models, XTENT and Gravity, and the two di erent ways to characterize the importance of settlements; based on population and the settlement's lifetime or solely on available archaeological data provided as input to the model.

Input data
As already noted, an agent may migrate only to a cell where known archaeological habitation sites existed in the specific time step, based on the archaeological survey conducted in the specific geographic area. The specific geographical locations of known archaeological habitation sites locations (and types) are provided as input spatial data file (shapefile) and used by the model during the simulation to map known site locations to available migrate locations for the household agents.

Submodels
Besides the settlements trading submodel, two other submodels are included in our ABM: the social organization of household agents within a settlement, in particular the self-organization social behaviour, and the natural disaster submodel. Detailed descriptions for the last two submodels are provided in previous work and versions of our ABM. Please refer to Section for the trading submodel and to Chliaoutakis & Chalkiadakis ( ) for the self-organization algorithm and Chliaoutakis et al. ( ) for details on the volcanic eruption submodel.