Toy Story: Homophily, Transmission and the Use of Simple Simulation Models for Assessing Variability in the Archaeological Record

The interpretation of spatial and temporal patterns in the archaeological record remains a long-standing issue in the discipline. Amongst many methods and interpretations, modelling of ‘biased transmission’ has proved a successful strategy to tackle this problem. Here, we investigate a type of biased transmission, homophily, that is the tendency of individuals to associate and bond with similar others. In contrast to other social sciences, homophily remains underused in archaeology. In order to fill this gap, we develop six distinct variants of a well-established modelling framework borrowed from social science, Axelrod’s Cultural Dissemination Model. These so-called toy models are abstract models used for theory-building and aim at exploring the interplay between homophily and various factors (e.g. addition of spatial features such as mountains and coastlines, diffusion of innovations and population spread). The relevance and implications of each ‘toy model’ for archaeological reasoning are then discussed.


Introduction
David Clarke famously described archaeological knowledge as 'a sparse suspension of information particles of varying size, not evenly randomly distributed in archaeological and Richerson (1985: 135) recognised three main classes: direct bias, when 'one cultural variant is simply more attractive than others'; indirect bias, when an individual perceives and selects a more attractive cultural variant by reference to another one; and conformist (or frequency-biased) transmission, when an individual picks a cultural variant based upon its higher frequency.
It is noteworthy that, whilst the use of random vs. biased transmission models in archaeology is deeply rooted in evolutionary thinking, methodologically comparable approaches were developed from an independent theoretical point of view in mathematical sociology as early as the 1950s (Skvoretz et al. 2004). For instance, random and biased net theory proceeds through the creation of nodes tied together through either random or biased processes, whose properties and effects upon the network structure are then analysed. One such bias is homophily, defined as 'the principle that a contact between similar people occurs at a higher rate than among dissimilar people' (McPherson et al. 2001: 416). In homophily, the intensity and likelihood of successful interaction is directly linked to a measure or perception of similarity between individuals. In a way, it is related, but not similar, to conformist transmission as individuals exhibit a tendency to model their behaviour upon a particular group characterised by a higher frequency of specific cultural traits. However, in homophily, these specific cultural traits can be confined to a given sub-group, and thus do not need be the most frequent and widespread across the entire population. Furthermore, the preference for a particular cultural variant in homophily is not simply based on its sole frequency, but on pre-existing similarity between individuals. Homophily and its structuring role in either simulated or real, observed networks are subject to a wide-ranging literature exploring various factors involved in this general process (see review in McPherson et al. 2001; from a mathematical point of view, see Fararo and Skvoretz 1987; for an evolutionary perspective: e.g. Fu et al. 2012, Haun and Over 2013, Creanza and Feldman 2014, Massen and Koski 2014. As any exhaustive review of the topic is beyond the remits of this contribution, we focus here on two properties of homophily. Firstly, it is important to recognise that, if homophily characterises the functioning of a particular social network, it does not necessarily apply to the entire corresponding society (e.g. Zhou 2011). For instance, Maoz's analysis of international political alliance and trade networks through the nineteenth and twentieth centuries AD shows that only the former are affected and shaped by homophily (Maoz 2012). Secondly, previous research has stressed that the impact of homophily is strongly related to the existing internal and spatial structure of the studied population (e.g. Kleinbaum et al. 2013;Ramazi et al. 2016). In this sense, it is worth introducing a distinction between, on the one hand, choice homophily where the interaction is driven by meaningful decisions taken by individuals and, on the other hand, induced homophily where interaction arises as a consequence of restricted cultural variation linked to population structure (e.g. McPherson and Smit-Lovin 1987;Kossinets and Watts 2009). Under such circumstances, homophilic interaction reinforces both existing intra-group cultural similarities and inter-group differences (Axelrod 1997; see also Centola et al. 2007).
By contrast, homophily remains relatively overlooked in archaeology, especially in the more modelling and quantitative-oriented literature. A recent paper by Shennan and colleagues provides a noticeable exception (Shennan et al. 2014). Drawing upon their statistical analysis of datasets of European Neolithic ornaments and pottery decorative styles, they demonstrate that the transmission process at play in pottery decoration was biased by a certain degree of homophily, although the same does not apply to ornaments, thus echoing known general properties of homophily (see above) and pointing to the presence of distinct historical trajectories within archaeological assemblages (Shennan et al. 2014: 108). These findings echo a large body of evidence drawn from a distinct archaeological sub-discipline with a long pedigree of interest in cultural transmission, namely ethnoarchaeology. Numerous studies have documented that distinct material traits are transmitted along different channels, thus implying varied forms and durations of social interaction (e.g. Pétrequin 1993;Gosselain 2000). In particular, techniques, such as pottery-making or flint knapping, are closely related to homophily, as they 'are socially learned and culturally transmitted […]. Their mastery requires learning through long lasting contacts, generally with socially close relatives' (Roux et al. 2017: 320).
The present contribution aims to fill this gap in archaeological method and theory by exploring homophily using a well-established modelling framework (Axelrod's Cultural Dissemination Model: Axelrod 1997;hereafter CDM). After reviewing the fundamentals of the CDM and its few archaeological uses, we develop a series of 'toy models' to explore the interplay of homophily and specific model variants, namely geographical space, diffusion of innovations and population structure. By labelling our simulations as toy models, we do not endeavour to denigrate their quality or relevance. On the contrary, following established terminology in numerous scientific fields, we stress the fact that these toy models are willingly and explicitly simple and abstract in design and aimed at theory building rather than grand narratives (see Alicea and Gordon 2014). The implications of toy models and self-advocated methodological simplicity for archaeological reasoning are briefly addressed in the final discussion and conclusion.

Axelrod's Cultural Dissemination Model
The CDM was elaborated 20 years ago by the political scientist Robert Axelrod to explore notions of convergence (Axelrod 1997). Axelrod's research question originated from the apparent paradox that, despite a linear relationship between the levels of interaction and the similarity shared by various individuals, global homogeneity never reaches completion. Thus, even if homophily is the main drive underlying interaction, the final population does not become 100% homogeneous. To solve this riddle, Axelrod developed a simple agent-based model.
The model comprises a population of agents arranged on a regular lattice and possessing several 'cultural' attributes. Agents exchange these traits with their neighbours through a process of homophilic interaction. To this purpose, each agent possesses a fixed number of features, which can be conceptualised as distinct cultural practices (e.g. language, technology). Each feature is then assigned a unique trait, taken from a list of possible variants. In a simulation based on five features (labelled A to E), for which there are 10 possible traits (labelled 1 to 10), an individual agent would thus be described as a list of, for instance, [A2, B6, C1, D5, E9]. The numbers of features and traits are set at the beginning of the simulation and do not change during each run. It must be noted that all traits and features are neutral and independent; the outcome of the interaction is independent of the actual content of each feature and only related to the overall level of similarity between agents. The model operates within a finite grid, each node corresponding to an agent. Each agent can only interact with its immediate neighbours. At the start of the simulation, this grid is completely populated with agents, who are randomly assigned traits for each feature. The simulation then proceeds through a repeated sequence of steps: 1. Randomly select one agent-known as source-on the grid 2. Randomly select one of its neighbours-known as target-to interact with 3. Calculate their cultural similarity as the number of features for which source and target have the same trait. With probability equal to this cultural similarity, source copies from target one feature, randomly selected from those features that source and target do not already share.
The simulation is initiated with each agent having a uniformly randomly selected set of traits, such that the initial population is very heterogeneous. Under these conditions, neighbouring agents have few traits in common, and their ability to interact together is relatively limited. As the simulation proceeds and an increasing number of interactions happens, the overall level of trait diversity drops, leading to the formation of 'regions', blocks of neighbouring agents sharing all traits, even though agents belonging to different regions are still able to interact together. The number of regions decays at a non-monotonic rate through the simulation to reach eventually a state of equilibrium when the grid is populated by a number of regions, whose respective agents are entirely different such that they can no longer interact together. The CDM demonstrates that interaction through homophily leads to both local convergence (i.e. the creation of homogeneous regions) and global divergence (as the regions eventually cannot interact together). In Axelrod's own words: 'The [CDM] shows how homogeneous cultural regions can arise without any intrinsic relationship between the separate dimensions that become correlated' (Axelrod 1997: 219). Axelrod's initial results demonstrate that the eventual number of regions at equilibrium depends on several parameters, most noticeably the initial number of features and traits (Axelrod 1997). A higher number of features lead to fewer, larger final regions, as this increases the probability of interaction. By contrast, a higher number of traits are paralleled by more numerous, smaller regions at equilibrium, as the larger range of existing traits reduces the probability of successful interaction. Two further elements affect the equilibrium state. Firstly, greater grid size increases the number of interactions needed to reach equilibrium, and these added interactions give more opportunities for regions to be assimilated, leading to reduced fragmentation of the final grid. Secondly, increasing the number of target source can interact with leads to a smaller number of regions at equilibrium, as the extended range of possible interactions provides more scope for convergence. It must however be noted that neither the precise number of resulting regions, nor their location, nor shape can be predicted in the standard CDM.
In the classical CDM, the number of features and traits remains constant during the simulation, a condition considered as highly unrealistic by Axelrod himself (Axelrod 1997). Klemm and colleagues demonstrated that the introduction of drift disrupts the typical progression of the CDM (Klemm et al. 2003). Drift is defined here as the given probability that, at each step, an agent will change a single trait by reverting to a value randomly selected from the existing repertoire. This probability is evenly distributed across the entire grid and independent of interaction. Klemm and colleagues' simulations demonstrate that such added drift introduces 'noise' in the CDM, which undergoes periods of stability and instability. At a certain level of drift, the CDM is disrupted and the system never reaches equilibrium (Klemm et al. 2003).
The CDM is a simple, robust and very successful model, as evidenced by its impressive legacy in social sciences and physics (Castellano et al. 2009). Perhaps in line with archaeologists' suggested relative lack of engagement with homophily and spatially explicit modelling, the CDM has only attracted limited attention by archaeologists. Two main exceptions can be found a few pages apart in a recent edited volume on cultural transmission during the Palaeolithic period. Madsen and Lipo's adaptation of the CDM lies in the aforementioned tradition of identifying forms of biased cultural transmission and aims at testing whether or not structured learning contributes to the formation of richer cultural repertoires observable in the archaeological record (Madsen and Lipo 2015). In this perspective, Madsen and Lipo introduced a hierarchy within traits, so that the acquisition of particular traits becomes a prerequisite for the successful transmission of other ones. Their results suggest a strong correlation between the diversity of cultural repertoires and high fidelity in teaching. Kovačević and colleagues modify the CDM to test the suggested identification of spatiotemporal patterns in European Upper Palaeolithic ornaments as ethno-linguistic groups (Kovačević et al. 2015). The changes made to the CDM are extensive as, for instance, agents are not set on a fixed location on a grid, but rather mobile using a random-walking algorithm. The interaction process is also altered. Agents are either conflicting, one of them-randomly chosen-imposing its entire set of traits to the next one, or sharing, and then pooling at random their traits. In one version of the model, the choice between conflict and sharing is randomly decided; in a second version, this choice is based on pre-existing levels of similarity, i.e. homophily. Lastly, both drift and mutation are incorporated in the model. Overall, simulations exploring a range of parameters show that the spatially structured distributions of archaeological ornaments can arise under both versions of the model (that is either random choice, or similarity-dependent), suggesting that these patterns could, but do not have to, be linked to past ethno-linguistic groups.
Beyond these specific examples, there are three reasons why the original CDM is of potential interest for archaeological practice. Firstly, the CDM offers a simple modelling framework entirely dependent upon homophily and spatial structure. Secondly, the CDM provides a robust, if serendipitous, analogy with archaeological assemblages and/ or 'cultures', feature and traits corresponding to any category of material culture and its variants (e.g. pottery and associated types). In contrast to many existing archaeological models focusing on a single type of evidence, the CDM provides a conceptual framework for exploring variability at assemblage level, and a test to David Clarke's 'assumption that the best model for archaeological entities is a polythetic model of some kind' (Clarke 1968: 37). Thirdly, the CDM is a spatially explicit model. This is particularly relevant given the aforementioned sensibility of homophily to spatial structure, although, admittedly, the CDM setting remains more abstract than the landscapes archaeologists are more familiar with.

Model Description
We implement six extensions of the CDM. Each of these toy models explores the interplay of homophily with a given parameter, and how this interplay can lead to a variety of patterns potentially observable in the archaeological record. The first two toy models ('mountains', 'coastlines') consist of the inclusion of spatial features in the homogeneous grid of the CDM. The following three models ('flow', 'new trait', 'new feature') explore the role of homophily in the diffusion of innovations. The last model ('wave') follows upon the last category by simulating the spread of a population and investigating the resulting effects on cultural variation.
All simulations described here work within a fixed grid, although the new variants create various discontinuities within this otherwise static space. Furthermore, whilst the CDM originally uses a square lattice, we replaced it with a hexagonal one as this reduces the isotropy of the grid, that is, the methodological artefacts that would occur due to spreading along a diagonal in a square grid (Randall et al. 2002). Sites are considered neighbours if they are immediately adjacent on this grid, with the majority of sites having six neighbours, except for agents located at the edge or in the corners which can have anywhere between two and five neighbours depending on their exact location.
All simulations were carried out in Matlab® version 2015b, utilising the Statistics and Machine Learning, Parallel Computing, and Distributed Computing toolboxes. Where summary statistics are used (e.g. arithmetic means and maxima in Fig. 4), these are calculated from 100 replicate simulations. Most simulations are run on a 21 × 21 hexagonal grid, with the exception of the wave toy model, which was run on a 12 × 48 grid. The original code for each toy model is available at the following online repository: https://github.com/NelisDrost/ToyModels.

Mountains
The first model is characterised by the addition of so-called mountains, which are nodes on the grid which are uninhabitable, devoid of any features and traits and with which interaction is therefore impossible. As a result, any agents located next to such barriers have less neighbours to interact with.
This toy model echoes previous work by Parisi et al. (2003), although with several key differences. Parisi and colleagues distributed agents and mountains in alternating parallel rows set in a square grid, thus reducing the interaction range for most agents to two neighbours. By contrast, our simulations explore a range of density and location of mountains (10%, 20%, and 40% of the grid size), which are uniformly randomly distributed. Our work presents other fundamental differences with Parisi and colleagues' contribution. Except for using a hexagonal grid, our simulation design follows the initial parameters set by Axelrod's 1997 paper (Axelrod 1997). On the contrary, Parisi and colleagues made several alterations to the conditions of the CDM. Firstly, features are binary and thus can only take two values (0 or 1). Secondly, they replace the CDM usual pairwise interaction by an 'assimilation rule' whereby agents change all of their features to match the more frequent values observed amongst their neighbours, thus replacing homophily by conformist (frequency-based) bias (Parisi et al. 2003: 165-6). Thirdly, their model also includes drift (see above). These numerous divergences not only guarantee the originality of our work but also preclude in-depth comparison of our results with those of Parisi and colleagues (see discussion).

Coastlines
The second toy model follows the same logic by introducing a spatial discontinuity within the grid. Several nodes, considered to be 'sea', are not allocated any feature or trait and cannot be interacted with directly. In contrast to the previous toy model, agents gain additional neighbours via 'sea routes' and are able to interact with both their immediate neighbours and cells which are less than three links away via sea nodes. As for the previous toy model, we explored a range of sea node densities (10, 20, 40%), uniformly randomly distributed on the grid.

Flow
Our third toy model introduces another spatial parameter to the CDM, in this case by changing the conditions for interaction between agents. Source is randomly chosen as per the usual CDM conditions, but there is a higher probability that it will preferentially choose target based on its direction from the source. In our example, agents preferentially interact with neighbours to their south-east. This preference is set evenly across the grid. The interaction process otherwise remains identical to the traditional CDM setup.

Innovation as New Trait
This toy model considers, under two distinct spatial variants, the effect of innovation upon the CDM by adding a novel trait to the existing repertoire. As for all other traits, the novel trait is neutral and does not increase the fitness of the agents in any way. Previous work by Tilles and Fontanari (2015) also modelled diffusion of innovation within a CDM framework, but our work differs in two ways. These authors restricted the presence of the novel trait to a single agent from the start of the simulation onwards. During the simulation, this novel trait could be acquired by other agents as per the normal homophily rules of the CDM, but it could not be lost at any point by the original agent. By contrast, we introduce the novel trait later during the simulation and amongst a subset of neighbouring nodes, by analogy to a process of locally distributed innovation. This decision is to imitate the structure of archaeological data as, in most cases, it proves impossible-if not futile-to track an innovation to a single point in time and space. More often, archaeological data rather allow us to highlight particular regions and periods where innovation occurred, as for instance various areas of plant and animal domestication across the globe (e.g. Fuller et al. 2014;Larson et al. 2014). The innovation time is changed across simulations using a logarithmic scale (innovation time set at steps 1, 16, 256, 4096,~600.000, 1.000.000). Time dictates the spatial scale of the innovation as the novel trait is introduced in a number of agents corresponding to the average size of regions at that moment. Two options are considered. In the first one, the novel trait is introduced in neighbouring agents belonging to an existing region (i.e. already sharing all traits across all features). In the second option, the novel trait is added amongst contiguous agents distributed in a random shape bridging different regions (hereafter referred to as 'block').

Innovation as New Feature
The rules dictating this variant are identical to the previous model but applied to a feature. The new feature is introduced at different moments during the simulation, either in an existing region or in a 'block' (see above). For the sake of simplicity, this new feature is associated with a single, neutral trait. This new feature is transmitted under a unique set of conditions, based on the general rules of the CDM: (1 and 2) a random source and target agent are selected as per the usual CDM rules and (3) cultural similarity is calculated only for those features that agents share: if one agent has only the five basic features, and the other agent has a sixth novel feature, only the first five features are compared. As before, interaction happens with a probability equal to their cultural similarity. If the source and target have the same set of features, or the target does not have the novel feature, a trait is copied at random as usual. If the target agent has the novel feature and the source does not, the source will randomly select from the novel feature and any features which it does not already share, each having equal probability.

Wave
The last toy model simulates the effect of an expanding population upon its culture. This requires modifications of both the spatial properties of the CDM, here affecting the original distribution of the agents, and the homophilic interaction process. At the start of the simulation, only a contiguous portion of the grid is populated with the usual random distribution of traits, and the rest of the grid is left empty. Interaction between agents occurs as per usual, with a single caveat: if either source or target is empty, no interaction occurs. In order to simulate a population expansion, a supplementary step is introduced. At a rate of D times per interaction steps, a source and a target agent are selected according to the normal rules. If the source is occupied and the target unoccupied, then source propagates in the empty target through cloning. It is worth pointing out that Parisi et al. (2003) independently attempted a relatively similar modification of the CDM, though in a distinct way. In addition to the already noted different parameters (i.e. replacement of homophily by conformist bias and introduction of drift), our respective simulations differ in two further ways: Parisi and colleagues only populate a single cell in an otherwise empty grid at the beginning of the simulation; in their setup, expansion can happen at each step, into all unoccupied cells neighbouring occupied cells.

Results
Although the toy models developed here all adopt the same core rules inherited from the original CDM, they each explore markedly distinct processes. As a consequence, rather than developing and using a standardised descriptive tool (e.g. region size and number are often privileged in the classic CDM literature), the results of each model are reported using distinct techniques and measures, so as to tackle both the distinct behaviours of each toy, and to address the intuition that first motivated each approach. While this approach precludes ready comparison between the toy models, our purpose is to better identify and understand the changes in behaviour resulting from individual models, and this is best achieved by adapting our analytical approach to each case. Figures 1, 2, 3, and 5 show maps similar to those found in Axelrod (1997), where the thickness of the line separating agents indicates the level of cultural dissimilarity between them (thicker is more dissimilar), and no line indicates complete similarity and membership of the same region. Figures 1 and 2 additionally show the presence of 'mountains' (black cells) and 'seas' (blue cells), as well as heatmaps of interaction rate, which count the number of successful interactions each agent undergoes (how many times the agent succeeds in copying a trait from a neighbour). Figure 1 shows the outcome of three individual simulations, corresponding to distinct densities and locations of mountains. Although these maps correspond to selected individual runs, they are representative of the results consistently observed across all simulations. As illustrated by the heatmaps on the right panel, the main impact of adding mountains to the CDM is a localised reduction in cultural instability-the number of times agents change traits over the course of a simulation. The amplitude of this effect changes according to the proximity and density of mountains. For instance, enclosed or isolated areas undergo far less changes before reaching a stable state. Conversely, narrow bottlenecks between two or more areas can undergo a very high rate of change, as the regions of stable culture that form in the neighbouring areas will have difficulty pushing through this constricted space.

Mountains
Another behaviour highlighted by this toy model is the importance of anchoring points to the formation of fixed boundaries between regions. To become stable, a region must have a fixed border with all its neighbours. Such border must either form a continuous loop or terminate at both ends at an anchor point, either against the edge of the grid, another stable region or a mountain. Adding mountains to the grid increases the availability of such anchoring points, allowing more stable regions to form and coexist.

Coastlines
As illustrated by Fig. 2, adding coastlines to the CDM leads to an increase in the number of successful interactions for agents located in the vicinity. This effect is due to the increased mobility of traits provided by sea travel and the resulting increase in diversity of traits to which each agent located near a 'coastline' is exposed. These simulations also tend towards a less diverse equilibrium, both because the grid is smaller (i.e. it presents less occupiable nodes) and because the added sea nodes reduce physical distance between agents, thus promoting cultural merging. As for the mountains toy model, whilst these effects are recurrent across all simulations, their implications for the exact size and shape of the resulting regions are not predictable.

Flow
As for the previous toy models, the results of this variation are best illustrated by a map of a single selection (Fig. 3). The introduction of directionally biased transmission leads to an alignment of fixed borders with the flow direction, as boundaries which run orthogonal to it are rapidly encroached upon. It is noteworthy that, if the overall shape of the regions can thus be predicted to some extent, the same does not apply to their exact size and/or location on the grid.

Innovation as New Trait
The results of this toy model are presented on Fig. 4, where graphs summarise data for a hundred runs each for both the introduction of the novel trait in a region or a 'block'. Under both scenarios and for all innovation times, very low mean values indicate that the novel trait does not spread extensively and disappears during the vast majority of simulations. This process however presents some variation. In a few instances when the innovation time is set after equilibrium has been reached (step 1.000.000; light blue line), the novel trait persists through the entire simulation, reaching ubiquity or nearubiquity. Another exception is noticeable: the much higher mean value for introduction at step 256 (yellow line) in a random block is skewed by a single case where the novel trait became associated with a region which became fixed at equilibrium.
For both innovation across regions and blocks, an increase in the persistence of the novel trait and the number of agents possessing it is linked to its gradually later introduction, although not in a linear way. Indeed, for the introduction of a novel trait at the beginning of the interaction in a given region (step 1, red line), this novel trait persists, at a low frequency, for a longer duration than novel traits introduced at later stage (steps 16, 256 and 4096, respectively, dark blue, yellow and purple lines).

Innovation as New Feature
The results of this toy model are not reported graphically. In the vast majority of cases, the novel feature quickly becomes ubiquitous across the entire grid. No difference in this behaviour between initial spatial distributions (region or block) could be discerned. The rate of expansion of the novel feature seems to be linked to the time of introduction, with expansion being slowest early on in the simulation, when successful interactions between agents are generally rarer due to the initial state of heterogeneity.
Wave Figure 5a compares the behaviour of a non-modified CDM to the behaviour of the 'wave' toy model, with a rate of expansion D = 1/5. As the original population expands, large homogeneous regions are formed at the expansion front. These large regions eventually break down, leading to a phase of increased cultural diversity. Over time, this diversity eventually decays, as per the usual behaviour of the CDM, until the system reaches equilibrium. Figure 5b shows this trajectory by plotting the diversity of different segments of the grid against time. Whilst the CDM literature generally uses measures of region size or count to describe the model, here, we use Shannon's diversity index, or Shannon's H, to track potential changes in the diversity of the cultural repertoire across time and space. This index, routinely used for instance in numerical ecology (Hill 1973;Borcard et al. 2011), is calculated for sites grouped in eight vertical slices across the grid. As the front progresses, its gradual loss of diversity is indicated by a consistent drop in the Shannon index values for the first total inhabitation of the corresponding new slice. Local diversity then rises to reach a local peak, corresponding to the gradual breaking down of the large regions associated with the expansion front. After this local peak, the toy model reverses back to the usual CDM behaviour, with a decay of the diversity across the entire grid until equilibrium is reached.

Discussion
The six toy models developed here present varying implications for both our general understanding of the CDM and the assessment of homophily in archaeological reasoning. The following discussion focuses on the latter.

Homophily and Space: Mountains, Coastlines and 'Flow'
As previously mentioned, homophily is known to be sensitive to population and spatial structure. In this perspective, our first two toy models introduce a further spatially explicit component in the CDM through the addition of either mountains or coastlines. These toy models work in opposite ways by respectively diminishing or expanding the quantity of neighbouring agents available for interaction. The effect of these variants is consistent as the number of time some agents change traits either decrease or increase, based upon the distance between these agents and mountains or coastlines. This effect remains spatially localised and does not exhibit any systematic implication regarding the number, position or shape of the regions. For instance, as coastlines favour interaction, they also prevent the formation of stable regions in their vicinity as agents have more opportunities to interact with other agents possessing different traits. The best way to describe this process is thus not through a measure of long-term similarity (i.e. stable regions) but through a measure of intensity of flow (e.g. heatmaps provided on Fig. 1). The implications of these two toy models for archaeological reasoning are most salient when considering the recent literature on networks. Several studies have shown the importance of the location and number of neighbouring nodes possessed by individuals, households or archaeological sites in explaining the structure and functioning of social and exchange networks (from an archaeological point of view, see Collar et al. 2015). Such nodes are generally identified and quantified through a measure of similarity of the archaeological record informed by more 'traditional' methods such as typology. The results of our first two toy models however stress the need to distinguish the technique used to describe the network (i.e. measure of similarity) and the eventual shape of the corresponding network. Under homophily in a heterogeneous physical space, sustained interaction between agents indeed does not translate, over time, into a culturally homogeneous network.
By contrast, if the intensity of interactions is directionally biased, our third toy model (flow) indicates that its direction impacts upon the shape of cultures. This recurrent spatial effect is self-reinforcing as, by skewing the probability of interaction towards one particular neighbour, one increases the overall similarity and thus the probability of future successful interaction. In a way, the simplicity of this toy model echoes an ideal archaeological situation where one area is, for whatever reason, perceived as a constant source of inspiration by its neighbours. This being said, as the effect of the preferred interaction direction is applied globally across the grid, the behaviour of this third toy Fig. 1 The left panel shows, after 500 k steps, the existing regions, the thickness of the delineating lines being representative of the proportion of traits shared between neighbouring agents (no line = all traits identical; thickest line: all traits different). The right panel shows the number of times an agent has changed traits, visualised as a colour ramp (red corresponding to low values, yellow to higher values). Density of 'mountains' equivalent to 10% (top), 20% (middle) and 40% (bottom) of the grid size model is qualitatively different from, for instance, core-periphery models favoured in archaeology. In this case, the biased flow of cultural influence may appear linear at a localised scale, as in our toy model, but at a regional scale, it may be found to be radiating from one or more cores.
Homophily and Innovation ('Innovation as New Trait', 'Innovation as New Feature') The next two toy models investigate the interplay between homophily and innovation. When a new trait is introduced, its success is limited in time and space, as it rarely persists when the simulation reaches equilibrium. The underlying mechanism is simple. As the novel trait is introduced amongst a small number of agents compared to the size of the grid, there is a high probability that its frequency is low relative to the overall repertoire of traits for the corresponding feature. Under these conditions, the novel trait is outnumbered and unlikely to persist through repetitive successful interactions between agents, especially when the innovation occurs towards the beginning of the simulation. Since traits are distributed in a random way at the beginning of the simulation, the probability of having two neighbouring identical agents is almost null, and the average size of a region is equal to one. The success of the novel trait is therefore in constant jeopardy pending upon the direction of a successful interaction. If the 'innovative' agent is the source of interaction, there is a probability that the novel trait will be eradicated through acquisition from target; if the interaction happens in the opposite direction, the novel trait will spread. Under such pressure, the potential for extensive dispersal of the novel trait remains very low, as reflected by mean value lower than one.
Several simulations however show that, sometimes, the novel trait persists for a longer period, though not until equilibrium. The duration of this persistence is driven by differences in innovation time (i.e. the later a novel trait is introduced, the greater its initial abundance) and, to a lesser extent, by the spatial mode of introduction of the novel trait (i.e. region vs. block). If the novel trait is introduced in a pre-existing region, it does not affect the ability of the agents to interact together, as it is simply substituted for an existing Fig. 2 The left panel shows, after 500 k steps, maps of the 'coastlines' toy model with a density of coastlines equivalent 10% (top), 20% (middle) and 40% (bottom) of the grid size. The right panel presents the corresponding heatmaps Fig. 3 Left: map of 'flow' toy model, with interaction biased to be more likely with a target south-east of the source, after 132 k steps. Right: map of the same flow toy model after 196 k steps shared trait. Yet, it does negatively affect their ability to interact with neighbouring agents belonging to other regions as the probability of sharing the novel trait with them is null. If the novel trait is introduced in a block, it increases the level of similarity between these agents, and thus the probability of their successful interaction, as they may not have already shared identical traits for this feature. As per the first option, the introduction of this novel trait reduces the probability of successful interaction with other neighbouring nodes. Therefore, pending upon the spatial mode of introduction, the innovation event is either detrimental or favourable to the local conditions of interaction.
Given that the new trait is neutral, its relative success does not depend upon any fitness or adaptive capacity. It is rather related to population structure, either reflected by the topology of the entire grid or by the distribution of the innovative agents. In nearly all instances, the spatial diffusion and temporal persistence of the novel trait remain limited, as their overall low frequency makes them highly susceptible to Fig. 4 a Plot of simulation time (x-axis, using a logarithmic scale) against mean count of agents with novel trait, when introduced in a region. b Plot of simulation time against maximum count of agents with novel trait, when introduced in a region. c Plot of simulation time (x-axis, using a logarithmic scale) against mean count of agents with novel trait, when introduced in a block. d Plot of simulation time (x-axis, using a logarithmic scale) against maximum count of agents with novel trait, when introduced in a block 'extinction' with the framework of the CDM. This behaviour recalls in many respects one facet of the European Palaeolithic archaeological record, characterised by a multiplicity of regional innovations in lithic technology. These present a wide range of outcomes, from limited success to long-lasting legacy. Several scholars have suggested that these multiple trajectories are not simply related to the adaptive potential of these technological traits but to levels of connections of the corresponding local populations (e.g. Hopkinson 2011, Malinsky-Bulller 2016, an interpretation echoed by the results of our toy model. When introduced as a novel feature, the innovation diffuses across the entire grid, regardless of the time or spatial mode of its introduction. The reason for this behaviour is simple. Whilst our toy model includes a mechanism for acquiring a new feature, there is however no process to lose it so that any successful interaction involving the novel feature contributes to its eventual global success. If simulations all reach the same outcome, the actual process differs whether the novel feature is introduced in a region or in a block. In the first option, the novel feature does not change the level of homogeneity of the region, although the probability of further successful integrations is now lower, as it is now calculated upon six rather than five features. In the second option, the addition of a new feature increases homogeneity amongst the involved agents, but this comes at the limited expense of future successful interaction, now distributed amongst six, rather than five, features. As these changes only affect the agents with the novel feature, it is not surprising that we did not observe any significant difference between both processes. Yet, there are potential implications, best perceived when considering introduction of the novel feature after equilibrium. If the novel feature is introduced into a region, the latter has already reached a fixed state, and the novel feature will thus remain confined to it. On the contrary, if the novel feature is introduced across different regions, its occurrence allows for agents previously totally dissimilar to interact again, restarting the homophilic interaction process. This discussion is not purely theoretical. The simplicity of this toy model provides a necessary baseline for the evaluation of further parameters. Such new variables could include adding trait diversity to the novel feature in order to see how it affects the spread of the feature and the general behaviour of the CDM. The interaction process used here is purely homophilic and in line with the normal rules of the CDM. Changes to this mechanism could be considered through, for example, the introduction of similarity thresholds required for successful transmission of the novel feature. In this sense, this toy model has probably more potential to assess how homophily shapes the diffusion of innovation (see also Skvoretz 1987: 1204-5).

Homophily and Population Structure: Wave
The wave toy model is perhaps the most immediately comparable to 'real-life' archaeological examples. In this toy model, the expansion of the new population in an empty landscape is accompanied by a loss of cultural diversity, reflected spatially by the creation of few, large homogeneous regions behind the wave-front. This trajectory is the direct outcome of the addition to the normal CDM of a step of propagation of the agents through cloning. This process leads to a loss of diversity, as agents located at the front exchange the same small range of traits. After the front has stopped its expansion, the toy model enters a second phase, characterised by a rise in diversity of the then entirely inhabited areas. The diverse and numerous traits present in the zone of origins eventually progress across the entire grid through homophilic transmission, at the same time breaking upon the initial homogeneous regions.
Although it is not possible to make any direct correspondence between a single step in the model and any given measure of time (e.g. day, season, year, generation), this two-phase trajectory can however be observed in many archaeological situations. Here, we only briefly discuss its relevance for the spread of early farming across Europe (see also Drost and Vander Linden in prep.). For example, Colledge and colleagues have previously identified a gradual loss of diversity in crop packages during certain stages of the spread of early farming across Europe (Colledge et al. 2004(Colledge et al. , 2005. Existing competing interpretations, based on various kinds of mathematical modelling, focus on the initial loss of diversity and suggest that this process is either related to biased transmission or random spatial drift (Conolly et al. 2008;Pérez-Losada and Fort 2011). Our results partially concur with the first explanation, as the loss of diversity is linked in our toy model to cloning, which can be conceptualised as an extreme form of biased transmission. But, in contrast to previous work, our toy model also accounts for a subsequent rise in diversity, a process observed in several distinct European regions and for several categories of evidence (e.g. crop packages in Neolithic Iberia: Antolín et al. 2015; Adriatic Neolithic zooarchaeological assemblages: Gaastra and Vander Linden 2018). Lastly, it is also worth noting that these results are robust to a wide range of parameters, including varying rates of dispersal, levels of drift, and proportions of local populations (Drost and Vander Linden in prep.).

Conclusion
Homophily is a powerful force of social interaction, yet it remains relatively overlooked in archaeology. As demonstrated by an extensive modelling and sociological tradition, for all its apparent simplicity, homophily is a highly variable process, shaped by a multiplicity of factors, in particular spatial and population structure. In order to explore this variability, and its potential implications for archaeology, we developed six toy models within the framework of the well-established CDM. Whilst this original model is based on simple rules and a continuous abstract space, the common thread of our diverse toy models is the creation of spatial and/or temporal discontinuities within the CDM. Furthermore, apart from the wave toy model, we did not amend the homophilic interaction process that underlies its functioning.
It is worth reminding, as many have done before us (e.g. Gerbault et al. 2014;Romanowska 2015), that models are not depictions of reality, but formal thought experiments designed to test and explore specific research questions, a position explicitly acknowledged by our focus on toy models. At the risk of repeating methodological clichés, simplicity in modelling is paramount for-at least-two reasons. Firstly, the sequential addition of single modifications allows for the precise evaluation of their effects and building basic knowledge required for understanding their interplay with further factors. Secondly, simple models allow preventing them from becoming 'black boxes', giving the impression that their functioning is difficult or impossible to assess, a perception detrimental to the use, understanding and dissemination of this approach, especially for outsiders (Cegielski and Rogers 2016).
Such 'simplicity by design' is reflected by the structure of our individual toy models, each linked to a specific variant to the CDM. The relevance of such parsimonious approach is best exemplified by the 'new feature' model which, as such, does not provide a lot of information but constitutes a robust baseline for further experiments. All toy models otherwise highlight the ability of homophily in generating a range of patterns in the spatial distribution of neutral, uncorrelated features and traits. Whilst the wave toy model shows the role of demic expansion in shaping cultural diversity, both mountains and coastlines demonstrate how simple landscape features inform cultural distribution. It must be stressed that, in all cases, such patterns arise independently of any adaptive value attached to any trait, a point also echoed in our discussion of the 'innovation as new trait' model.
These toy models present archaeological potential because they contribute to theory building by offering possible alternative hypotheses to account for known archaeological situations, and because they all raise questions regarding the ways we interpret the variability of the archaeological record. There is a long-standing tradition in the discipline to reduce this variability to measures of similarity and corresponding spatial and temporal patterns, as epitomised by concept of archaeological cultures. Our results suggest that there is as much, if not more, information to be gained in frequency (i.e. number of times agents change traits) and diversity (i.e. variation in the cultural repertoire). In this perspective, an avenue worth exploring lies in the development of measures appropriate for both simulation and real-life data, although complications are likely to arise because of the different resolutions of such datasets (i.e. controlled simulation data vs. biased archaeological data; see, from a CDM perspective, Kovačević et al. 2015).
There is no doubt that the archaeological variability cannot be reduced to the joint effects of homophily and various spatial parameters. Yet, a simple, incremental method as advocated here is a more productive long-term strategy than a priori assumptions of complexity and holism separated from any formal analyses. In this perspective, focusing on homophily, or any other factor for that matter, should be seen as an addition to an already large array of interpretations, as well as an invitation to consider time and again new ways to conceptualise and describe archaeological variability.