Assortative pairing alone can lead to a structured biota in organisms with cultural transmission

Spatial separation is often included in models of ethnic divergence but it has also been realised that urban subcultures can, and frequently do, emerge in sympatry. Previous research tended to attribute this phenomenon to the human tendency to imitate self-similar individuals and actively differentiate oneself from individuals recognized as members of an outgroup. Application of such a model to non-human animals has been, however, viewed as problematic. We present a parsimonious model of subculture emergence where the algorithm of social learning does not require the assumption of an 'imitation threshold'. All it takes is a slight modification of Galton-Pearson's biometric model previously used to approximate cultural inheritance. The new model includes proportionality between the variance of inputs (cultural 'parents') and the variance of outputs (cultural 'offspring'). In this model, assortment alone can lead to the formation of distinct cohesive clusters of individuals (subcultures) with a low within-group and large between-group variability even in absence of a spatial separation or disruptive natural selection. Sympatric emergence of arbitrary behavioural varieties preceding ecological divergence may thus represent the norm, not the exception, in all cultural animals.


Introduction
There is a long-standing debate about the possibility of speciation that is sympatric, i.e., the emergence of multiple species from one species without geographic isolation [1][2][3][4][5][6][7][8] . Although the importance of the geographical aspect of speciation had lately been questioned 9,10 , spatial isolation is still believed to play a pivotal role in divergent evolution because it limits the gene flow between adjacent populations 11 . In this contribution, we investigate the implications of assortative mating within a model of nonparticulate multidimensional inheritance. The inheritance model can be used to approximate both an additive genetic transmission (Galton-Pearson model 12 ) or a cultural one (Parental Variability-Dependent Inheritance 13 ). We demonstrate that in a model where the effect of culture is dominant, assortment can even on its own, i.e., in absence of any disruptive natural selection, lead to the formation of distinct cohesive clusters of individuals with low within-group and large between-group variance. Parental Variability-Dependent Inheritance can explain the frequent sympatric emergence of subcultures and varieties characterised by limited gene flow (although the subcultures or varieties can diverge ecologically at a later point). This model might be particularly relevant for speciation in taxa with a high capacity for cultural evolution, such as the cetaceans 14 , songbirds 15,16 , or hominins 17 , which are known for their potential to undergo adaptive ecological radiation.

Darwin's mystery of mysteries
Originally, the famous designation 'mystery of mysteries' 18 was not reserved just for the formation of distinct species. It referred to the emergence of distinct groups of organisms with high between-group and low within-group variation 19 . Darwin stressed the absence of fundamental differences between varieties, subspecies, species, and higher taxa 20 . Differences between such categories were supposed to be the question of scale rather than quality. Current research tends to support this view 9,21 . It seems that reproductive isolation can emerge quite easily even in conditions characterised by an uninterrupted gene-flow 7 . Two distinct varieties may develop into two separate species via several broadly overlapping stages: (1) homogeneous gene flow between populations, (2) heterogeneous gene flow where alleles directly linked to diverging traits rarely cross the boundaries between populations, and finally (3) separate species that meet the condition of full reproductive isolation 22 . It should be noted, however, that resistance from the mainstream against this view is still quite considerable 11 .
Nevertheless, the subject of structured biota is no less challenging than the better-known speciesoriented perspective 19 . Why do, for example, distinct cultures, subcultures, or ethnic groups exist? This is difficult to answer because research in cultural evolution tends to focus on arbitrary signs or ornaments which do not necessarily have an adaptive function that would allow for an ecological explanation of the divergence 23 . Computer simulations of minimal counterfactual systems can shed light on this conundrum.

Nonparticulate inheritance models
Two models of nonparticulate inheritance have recently been employed in studies of cultural evolution 12,13 . These models do not rely on, in this context restrictive, beliefs about additional genetic variation: they can be viewed as general inheritance models such as those which preceded or existed in parallel with Mendelism. They have been overlooked due to evolutionary biologists' focus on genetic models which are believed to provide a good approximation of the simplest form of organismal inheritance. Both models are based on simple blending inheritance but replace the problematic assumption of the offspring ( ) being exactly in the arithmetic mean ( ( 1 , … , ) = 1 ∑ =1 ) of parental values ( 1 , 2 ) with the assumption of a random normal distribution of offspring (N( , ), with mean and variance ). Galton-Pearson (GP) model (Equation (1.1)) assumes a mutation term characterised by a constant standard deviation ( ) independent of the difference between parental values.
There, the offspring assumes a position on a vector connecting her parents. Her distance from the point in-between parental positions is normally distributed. The lower the proportional offspring variation ( ) relative to the constant offspring variation ( ), the closer is the model to a pure Galton-Pearson inheritance. A visual comparison of the models can be found in Figure 1.  2)). In a system with pure PVDI, the proportion of offspring between parental values is constant. In a system with GP inheritance, the proportion depends on distance between parental values.
Assortative social learning and mating seems almost omnipresent in human populations 25 and has been abundantly documented also in nonhuman animals 16,26 . It has been previously identified as a necessary but not sufficient condition of sympatric speciation 3 .

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. The difference between two individuals can be calculated as the distance between their positions in a trait space 27 . We assume no external natural selection acting on the population and there are only relative preferences for pairing between individuals. In our simulations, relative preferences between agents are determined by a homophily coefficient (ℎ) which expresses a preference for self-similarity. If ℎ is negative, we observe a negative assortment where dissimilar individuals are more likely to pair. When ℎ is positive, there is a positive assortment and similar individuals pair with a higher probability. For ℎ = 0, there is no assortment, and the pairing is completely random.
We can picture this as a cultural transmission or procreation in a multidimensional trait space. In each step, each agent selects an interaction partner ('role model') and modifies her position along the vector between her and her interaction partner. An agent will probably move towards the partner, sometimes overshooting her original position. With the same probability as overshooting her position, the agent moves away from the partner (in Figure 1, 1 can be read as the agent's original position, 2 as the selected role model's position, and offspring density as the agent's resulting position distribution after the interaction).

The formal model
In each simulation run, each agent in a population of size is represented by its position in a Euclidian -dimensional space (a 'trait space'): The difference between any two individuals in a trait space can be calculated as the distance between two points in -dimensional space: 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' , each agent selects a partner probabilistically. Relative preference of the focal individual for agent is given by where ℎ is a homophily coefficient, that is, a measure of preference for self-similarity. This wellbehaved function ensures that self-preference, i.e., preference for individuals whose distance from an agent is zero, serves as a referential value of 1, and other preferences follow proportionally depending on ℎ.
The probability that selects is given by If ℎ is negative, we observe heterophily, that is, negative assortative mating where dissimilar individuals pair with a higher probability. For positive ℎ, pairing follows homophily, a positive assortative mating, and it is similar individuals that pair with a higher probability. For ℎ = 0, there is no assortment, relative preference for all individuals is 1, and pairing is completely random. The absolute value of homophily |ℎ| corresponds to the strength of assortment. For ℎ = 1, agent who is twice closer to the focal agent than agent will be selected as a partner by approximately twice as often as (assuming ‖ − ‖ and ‖ − ‖ are significantly larger than 1). For ℎ = 2, will be selected approximately four times as often as . For ℎ = 3, will be selected approximately eight times as often as , etc. Similarly, for ℎ = −1, will be selected twice as often as etc.

mutatis mutandis.
Because each individual ( ) at time step selects one role model ( ), we obtain pairs of parents in all time steps. Each pair of parents creates one offspring, denoted for this purpose ( + 1), whose position in the trait space is given by probabilistic non-particulate inheritance combining Galton-Pearson and PVDI terms.
Agent's position in the following + 1 timestep, ( + 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 S4). The distance of ( + 1) from the arithmetic average of ( ) and ( ) is normally distributed.
The construction of agent configuration in + 1, ( + 1) from ( ) takes only a single time step, which means that all agents in our simulation alter their positions in synchrony. The 'inspirational' model, which is the focus of this manuscript, does not require any additional specifications.

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 = 10.
From the initial configuration, the partner selection + position inheritance algorithm was iterated for 200 steps in each simulation run (400-800 runs per parameter combination depending on the detail required given the variation in the number of distinct clusters).

Analysis and visualisation
The average distance between agents in the trait space was calculated to assess whether the population In the main manuscript, we abstained from any extensions, including normalisation, because our aim was to demonstrate the potential of a very simple model 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 S3.
Demonstrative examples were run using R 28 . The code for parallel runs was written in Python 29 using Numpy 30 infrastructure and delegated to processor cores using R package parallel 28 and reticulate 31 .

Results
We found that structured biota, manifested as distinct groups with high between-group and low within-group variance, emerges from positive assortment alone when PVDI prevails over GP, that is, when the constant standard deviation in an inheritance model that combines the GP with PVDI is small ( Figure 3B). The emergence of distinct subcultures was observed even when -in contrast to some previous studies 32 -we did not assume that agents coordinate their positions below or extremise their positions above certain threshold distance.
In a system with a considerable influence of constant offspring variation, assortment alone cannot lead to the emergence of a structured biota ( Figure 2). Instead, the agents form a single uniform cluster ( Figure 3A). This outcome is altered neither by homophily nor by population size.
Regardless of the constant offspring standard deviation , smaller populations face a reduction of effective trait-space dimensionality (Supplement S1, Figure S1). The effective dimensionality decreases when relative offspring variability is high. In smaller populations, unidimensional polarisation or variability loss is more likely, which conforms to the results of previous research 33  The inability of the Galton-Pearson system, which approximates the model of polygenic additive inheritance, to form distinct varieties in the presence of assortative pairing alone might be interpreted as providing support for ecological theories of sympatric divergence 7,9,11,18 . In absence of any cultural inheritance, natural and sexual selection must take place concurrently to generate separate ecotypes or species 4,8,18 .

Discussion
We argue that unfragmented species of sexually reproducing organisms with inheritance system that follows the PVDI (such as culture) are held together by the same force that holds together the morphological lineages of asexual organisms, that is, by stabilising natural selection. In such species, sympatric speciation should not be viewed as an exception but rather as a norm. Distinct variants that emerge spontaneously can easily settle into distinct ecological niches if there is a potential for divergent selection 14 . Such clusters are precursors of cultures, subcultures, guilds, alternative subsistence strategies, political factions, etc. 34 . We suggest that this process enables organisms with the capacity for social transmission to form and maintain distinct ecotypes defined by rituals, ethnic markers, and dialects even prior to the emergence of any genetic differences between these groups 35 .
Such organisms should be therefore capable of exploiting natural resources in ways which are out of the reach of organisms that rely solely on genetic inheritance. Sympatric fission, that is, cases where cultural divergence precedes any spatial dislocation of individuals, may play a vital role in the distribution of cultural variants across time and space 14,34 . Gradually, one may see an accumulation of genetic differences aiding further differentiation between groups 22 . Cultural assortment sparked by preference for self-similar individuals may be the key factor that limits gene flow between emerging subpopulations. The attested historical adaptive radiations of hominins 17 may be in part due to the vast cultural capacity of our lineage.           The same set of simulation runs as in Figure S11 was used to generate the image. PC1-PC3 stand for the first three Principal Components. See the explanation of the dimension reduction in a system with large and small in the caption of Figure S14.  For detailed description see the caption of Figure 3 in the main article. Simulation run parameters were chosen to match those of Figure 3 Equivalents of Figures S3-S10 are not included for brevity, but can be easily obtained using the code at https://osf.io/pvyhe/?view_only=a79f0d07847f45cab069a8b8f09f258b. Notice the clustering tendencies and dimension reduction in the system with large influence of GP inheritance and homophily. The polarization is due to inflation of the originally small differences. If the homophily is large and there is normalization after each simulated generation, distances change proportionally to their magnitude. Similarly as in a system with PVDI and without normalization (see Figure 3B), larger distances get larger, which may cause polarization.     For detailed description see the caption of Figure 3 in the main article. Simulation run parameters were chosen to match those of Figure 3 Equivalents of Figures S3-S10 are not included for brevity, but can be easily obtained using the code at https://osf.io/pvyhe/?view_only=a79f0d07847f45cab069a8b8f09f258b.