2.1. Overview

We employ an agent-based simulation to study the emergence of a structured biota under the assumption of assortative pairing. Individuals close to each other in a multidimensional trait space are considered similar. We assume no external natural selection acting on the population and only relative preferences for pairing between individuals. In our simulations, relative preferences between agents are determined by a homophily coefficient (h) expressing a preference for self-similarity. If h is negative, we observe a negative assortment where dissimilar individuals are more likely to pair. When h is positive, there is a positive assortment and similar individuals pair with a higher probability. For h = 0, there is no assortment, and the pairing is completely random.

In each step, each agent selects an interaction partner (‘role model’) and modifies her position along a vector between her and her interaction partner. The value $t_{p_1}$ in Figure 1, can be in this context interpreted as the agent's original position, $t_{p_2}$ as the selected role model's position, and offspring density as the agent's resulting position distribution after the interaction.

We use the continuous inheritance model, with PVDI and GP terms combined, since a pure PVDI model implies that if cultural parents share a trait value, the offspring become perfect learners who cannot make errors. The standard deviation in the combined model is controlled by two parameters, constant offspring variation η and coefficient of proportional variability ν, that are introduced in the description of continuous cultural inheritance models above.

We offer this narration for the combined model: agent preferences are centred around a mean of their cultural parents and their standard deviation is proportional to the standard deviation of the cultural parents: N(μ(t _p), ν ²σ ²(t _p)), where the variance of the distribution is just the standard deviation squared. Because an agent is not always able to match its preferences exactly, the resulting position is the sum of the preferred position and some random influence independent of parental traits: N(0, η ²). The additivity of variance (Fisher, Reference Fisher1918), N(μ(t _p), ν ²σ ²(t _p)) + N(0, η ²)  ≈  N(μ(t _p), ν ²σ ²(t _p) + η ²), allows us to draw just one random number per modification of cultural position and yet consider both sources of randomness in cultural inheritance simultaneously (see Supplement 1 for details). With the probability proportional to the pink area in Figure 1, the agent will end up between its original position and the role model's position. For cases where η = 0, this probability p _between can be calculated from ν using

$$p_{\mathrm{between}} = \displaystyle{2 \over \pi }\mathop \smallint \limits_0^{{1 \over {\nu \sqrt 2 }}} e^{{-}t^2}dt, \;$$

which is derived from the normal cumulative distribution function (see Supplement S2 for details). Because normal distribution is symmetric, the agent extremises its own position away from the role model with the same probability as it overshoots the role model's position along the vector that connects them, therefore

$$p_{\mathrm{extremise}} = p_{\mathrm{overshoot}} = \displaystyle{{1-p_{\mathrm{between}}} \over 2}.$$

For η > 0 these formulas would depend on the distance between the agents. However, for pure PVDI, we can calculate that e.g. p _between = 0.68 for ν = 1, and p _between = 0.38 for ν = 2, or for many regularly spaced values of ν we find that p _between = 0.5 when ν = 1.48 or that p _between = p _extremise = p _overshoot when ν = 2.32. (The values of ν were obtained using numerical methods and rounded to 2 decimal places. Analytical solution is not possible, because the Gaussian error function has no closed-form expression.) On a similar note, it is possible to run a series of simplified simulations and find out that for ν ≥ 2.10 two points that serve as role models to each other are more likely to end up further apart after the cultural exchange than they were before it.

While the example with skirt length (or canoe size, spear length, ratio of red and white cattle, age at marriage, amount of milk to put in coffee, etc.) represents a special unidimensional case of our model, we focus on cultural transmission or procreation in a multidimensional trait space. We do not need to suppose that in reality, the trait space that is used to represent differences and similarities between individuals has such a straightforward one-to-one correspondence with the design-space of an artefact (Mesoudi & O'Brien, Reference Mesoudi and O'Brien2008). The model bears a strong resemblance to the idea of a broad underlying culture space described in the work of Sperber (Reference Sperber1996). It represents a useful tool that generalises to discrete cultural variants if the probability (continuous) that certain behaviour is executed becomes the centre of the analysis instead of the behaviour's single occurrence (binary).

2.2. The formal model

In each simulation run, each agent A _i in a population P of size n is represented by its position in a Euclidian D-dimensional space (a ‘trait space’ or ‘culture space’):

(4)$$P = ( {A_1, \;A_2, \;\ldots , \;A_n} ) , \;\;\;A_i = ( {t_{i, 1}, \;t_{i, 2}, \;\ldots , \;t_{i, D}} ) \;\epsilon \;{\rm {\opf R}}^D.\;$$

The difference between any two individuals in the trait space can be calculated as the distance between two points in a D-dimensional space:

(5)$$\l\Vert A_i-A_j \r\Vert = \sqrt {\mathop \sum \limits_{d = 1}^D {( {t_{i, d}-t_{\,j, d}} ) }^2} .\;$$

In this trait space, we are simulating a process of non-particulate position inheritance. There is no external natural selection operating on the population and only relative preferences for pairing between individuals. In each time step or ‘generation’ g, each agent selects a partner probabilistically. The relative preference of the focal individual A _i for agent A _j is given by

(6)$${\rm rel.pref}( {A_i, \;A_j} ) = \left({\displaystyle{1 \over {1 + \l\Vert A_i-A_j \r\Vert}}} \right)^h, \;$$

where h is a homophily coefficient; a measure of preference for self-similarity. This well-behaved function ensures that self-preference, i.e. a preference for individuals whose distance from an agent is zero, serves as a referential value of 1, and other preferences follow proportionally depending on h.

Relative preference of agent A _i for agent A _j divided by the sum of all relative preferences of A _i defines probability

(7)$${\rm prob}( {A_i, \;A_j} ) = \displaystyle{{\mathrm{rel.pref} ( {A_i, \;A_j} ) } \over {\mathop \sum \nolimits_{k = 1}^n \mathrm{rel.pref} ( {A_i, \;A_k} ) }}\;$$

that A _i selects A _j as role model. Role models are sampled with replacement; one agent can therefore serve as the selected role model to more than one agent.

If h is negative, we observe heterophily, that is, negative assortative “mating” where dissimilar individuals pair with a higher probability. For a positive h, pairing follows homophily, positive assortative “mating”, where similar individuals pair with a higher probability. For h = 0, there is no assortment, relative preference for all individuals is 1, and pairing is completely random. The absolute value of homophily |h| corresponds to the strength of assortment. For h = 1, agent A _j who is twice as close to the focal agent A _i as agent A _k, will be selected as a partner by A _i approximately twice as often as A _k (assuming A _i − A _j and A _i − A _k are significantly larger than 1). For h = 2, A _j will be selected approximately four times as often as A _k. For h = 3, A _j will be selected approximately eight times as often as A _k, etc. Similarly, for h = −1, A _k will be selected twice as often as A _j etc. mutatis mutandis.

Because each individual A _i(g) at time step g selects one role model A _j(g), we obtain n pairs of parents in all time steps. Each pair of parents creates one offspring, denoted for this purpose A _i(g + 1), whose position in the trait space is given by probabilistic non-particulate inheritance combining GP and PVDI terms (see Eqn 3), that is,

(8)$$A_i( {g + 1} ) = \mu ( {A_j( g ) , \;\;A_i( g ) } ) + \displaystyle{{A_j( g ) - A_i( g ) } \over {\l\Vert A_j( g ) - A_i( g ) \r\Vert }}{\rm N}\left({0, \;\eta^2 + \nu^2\displaystyle{{\l\Vert A_j( g ) -A_i{( g ) \r\Vert}^2} \over 4}} \right).$$

The agent's position in the following g + 1 timestep, A _i(g + 1), therefore always lies on a line connecting the parental points (a variant of the model with additional normal noise can be found in Supplement S7). The distance of A _i(g + 1) from the arithmetic average of A _i(g) and A _j(g) is normally distributed. This approach has the advantage of arriving at identical results under any rotation of coordinate basis.

The construction of agent configuration of a population P in g + 1, P(g + 1), from P(g) takes only a single time step, which means that all agents in our simulation alter their positions in synchrony.

A similar model that pairs individuals exclusively for each time step can also be constructed. With this modification, n/2 pairs of agents are formed in generation g, each creating two new agents, ensuring that the population size remains constant, whereby the subsequent time step g + 1 follows the same inheritance algorithms. This modification invites new interpretations: (1) a vertical cultural transfer from biological parents to biological offspring, meaning that one step is identified with a biological generation; or (2) exclusive interactions, such as discussions, conversations and exchanges of opinions, after which the positions of both interacting individuals change simultaneously. In this model, an additional specification of the algorithm must ensure within-pair exclusivity. The selection of interaction partners in step g is decided one agent at a time in a random order. If an agent is already selected as an interaction partner, she does not select a partner of her own and cannot be selected as and interaction partner by any other agent. In any case, both models – let us refer to them as ‘inspirational’ and ‘interactional’ respectively – are highly similar and lead to equivalent conclusions. The ‘inspirational’ model, which is the focus of this paper, does not require any additional specifications.

2.3. Computer simulation

The initial state of the computer simulation was set to a single cloud of normally distributed points with standard deviation σ = 100 along each trait value. Simulation results were obtained for D = 10.

Simulation parameters (n = 250, 500, 1000, η = 1, 10, 100, 0 ≤ h ≤ 4, 0 ≤ ν ≤ 3.0) were chosen so as to demonstrate important transitions between systems with different frequencies of divergence. From the initial configuration, the partner selection + position inheritance algorithm was iterated for 200 steps in each simulation run (1000 runs per parameter combination).

2.4. Analysis and visualisation

The average distance between agents in the trait space was calculated to assess whether trait space collapses into a single point (variability loss) or enters a feedback loop of ongoing expansion (variability explosion, see Supplement S4, Figure S3). The effective dimensionality of the trait space was quantified by the proportion of variance of agent positions explained by the first three principal components (PCs) of the D-dimensional trait space and by the number of PCs necessary for capturing 99% of all variance. These measures were highly correlated (cor = 0.9), so in the following we only worked with the first measure. We chose variance explained by the first three dimensions because it has an intuitive interpretation in the context of the resulting images: what proportion of overall variance in 10 dimensions is captured in the 3D scatterplot. If this proportion is for example 95%, we can conclude that despite the model running in 10 dimensions, the effective dimensionality of the population is lower because we would obtain an almost identical distance matrix capturing differences and similarities between agents if we used only the PC1–PC3 space. The number of distinct clusters was evaluated using the HDBSCAN algorithm (for a description of the algorithm, see Supplement S3), a hierarchical version of DBSCAN (Density-Based Spatial Clustering of Applications with Noise), which is an important part of the rich tradition of density-based clustering methods in biology (Edla & Jana, Reference Edla and Jana2012). Unlike DBSCAN, HDBSCAN sidesteps the problem of how close points should be to be considered a part of the same cluster.

The development of agent network in the trait space is visualised across time steps as a series of static images and as an animation capturing the dynamical process of clustering (see Supplementary animations). The 3D scatterplots work with the first three PCs (standardised to have mean = 0 and SD = 1 along each PC) and summarise the layout of points in a D-dimensional trait space. To minimise scatterplot rotation between adjacent frames, PC1, PC2 and PC3 are mapped to axes x, y and z through PC rearrangement and reversion, such that cor(x _g, x _g+1) + cor(y _g, y _g+1) + cor(z _g, z _g+1) is maximised. The points indicating agents’ relative positions in PC1–PC3 space are coloured according to their assignment to distinct clusters. The biggest subgroup from cluster Γ in time step g that belongs to a single cluster in time step g + 1 inherits the colour of Γ from the previous step. The measure of effective dimensionality (Figure S2), the number of distinct clusters (Figure 2), and the average distance between agents (Figure S3) are also visualised as variables dependent on time in a graphical summary of a single simulation run. The proportion of non-clustered noise is assesed at the end of each simulation run (Figure S4).

Figure 2. A graphical summary of the number of subcultures after 200 model generations. The points in the 10-dimensional culture space were normally distributed across all dimensions at the beginning of each simulation run, and 1000 simulation runs were executed for each parameter combination. Red crosses indicate the values of parameters used for single-run examples included in the main article. Blue crosses indicate examples available in the Supplementary Material.

The average values of measures for all simulation runs with the same parameter configuration are visualised as colour shades depending on parameters n, h, ν, and η.

A straightforward extension of the model allows us to focus on relative rather than absolute differences between individuals in processes that drive cultural acquisition. In such a model, we normalise the positions of agents so as to maintain the average distance between agents constant after each step. Thanks to this additional step, the model leads neither to variability loss (where ν and η are too low) nor to variability explosion (where ν is too high), which were both present in the original model. In effect, normalisation binds η to the overall population variance while keeping it independent of the two particular values selected as parental traits. This resembles Fisher's genetic model which reconciled biometric inheritance with Mendelism under the assumption of additive genetic variance (Fisher, Reference Fisher1918). In the main paper, we abstained from any extensions, including normalisation, because our aim was to demonstrate the potential of a very simple model, one relying only on bilateral relationships between interacting agents, to sustain a reasonable cultural variation. A more detailed exposition of the normalised model can be found in Supplement S6.

Demonstrative examples were run using R (RCoreTeam, 2019). The code for parallel runs was written in Python (Van Rossum & Drake, Reference Van Rossum and Drake2009) using Numpy (Van Der Walt, Colbert, & Varoquaux, Reference Van Der Walt, Colbert and Varoquaux2011) infrastructure and delegated to processor cores using R packages parallel (RCoreTeam, 2019) and reticulate (Allaire & Ushey, Reference Allaire and Ushey2020). Visualisations were created using the base R graphics.