Mathieu Forget
email: forget@bio.ens.psl.eu

Sandrine Adiba

Leonardo Gregory Brunnet

Silvia De Monte

Heterogeneous individual motility biases group composition in a model of aggregating cells

Aggregative life cycles are characterized by alternating phases of unicellular growth and multicellular development. Their multiple, independent evolutionary emergence suggests that they may have coopted pervasive properties of single-celled ancestors. Primitive multicellular aggregates, where coordination mechanisms were less efficient than in extant aggregative microbes, must have faced high levels of conflict between different co-aggregating populations.

Such conflicts within a multicellular body manifest in the differential reproductive output of cells of different types. Here, we study how heterogeneity in cell motility affects the aggregation process and creates a mismatch between the composition of the population and that of self-organized groups of active adhesive particles. We model cells as self-propelled particles and describe aggregation in a plane starting from a dispersed configuration. Inspired by the life cycle of aggregative model organisms such as Dictyostelium discoideum or Myxococcus xanthus, whose cells interact for a fixed duration before the onset of chimeric multicellular development, we study finite-time configurations for identical particles and in binary mixes. We show that co-aggregation results in

Introduction

Aggregative life cycles are set to the boundary between unicellular and multicellular organization. In aggregative life-cycles, initially sparse cells come together and form multicellular groups. Multicellular development following aggregation is essential for the formation and successive dispersal of stress-resistant dormant cells. Such cysts or spores seed subsequent generations, where cells undergo unicellular, vegetative growth [START_REF] Raper | Pseudoplasmodium formation and organization in Dictyostelium discoideum[END_REF]. Both in prokaryotes, e.g. Myxobacteria, and in eukaryotes, e.g. Dictyostelids, multicellular development is associated with cell differentiation and division of labor, analogous to organisms, like metazoans, whose body derives from clonal growth. These two orders, and in particular the model species Myxococcus Xantus and Dyctiostelium discoideum, have thus been extensively regarded as highly informative on the transition from unicellular to multicellular organization [START_REF] Strassmann | Evolution of cooperation and control of cheating in a social microbe[END_REF][START_REF] Forget | Social conflicts in Dictyostelium discoideum : a matter of scales[END_REF][START_REF] Vos | Social Conflict in Centimeter-and Global-Scale Populations of the Bacterium Myxococcus xanthus[END_REF][START_REF] Kraemer | Endemic social diversity within natural kin groups of a cooperative bacterium[END_REF].

A key feature of aggregative multicellular species is their motility, that allows cells to actively pursue their bacterial preys during the solitary phase of the life-cycle, and to engage in a series of collective behaviours, such as swarming [START_REF] Shi | Chemotaxis plays a role in the social behaviour of myxococcus xanthus[END_REF] and collective taxis [START_REF] Bonner | The orientation to light and the extermely sensitive orientation to temperature gradients in the slime mold dictyostelium discoideum[END_REF]. Individual cell movement moreover allows cells to rapidly coalesce into multicellular aggregates as a response to stressful environmental conditions, typically nutrient exhaustion. Attaining collective functions without delay has been argued to constitute an important adaptive advantage to aggregative versus clonal multicellularity [START_REF] Márquez-Zacarías | Why have aggregative multicellular organisms stayed simple?[END_REF].

As a side effect, it is almost inevitable that motility leads to mixing of heterogenous cells to a much larger extent than in colonies of sessile microbes or clonally growing multicellular organisms. Indeed, genetic chimerism occurs naturally in aggregative multicellular species [START_REF] Strassmann | Altruism and social cheating in the social amoeba Dictyostelium discoideum[END_REF]Fortunato et al., 2003b;[START_REF] Gilbert | High relatedness maintains multicellular cooperation in a social amoeba by controlling cheater mutants[END_REF][START_REF] Vos | Social Conflict in Centimeter-and Global-Scale Populations of the Bacterium Myxococcus xanthus[END_REF][START_REF] Kraemer | Endemic social diversity within natural kin groups of a cooperative bacterium[END_REF][START_REF] Sathe | Interspecies and intraspecies interactions in social amoebae[END_REF]. Within chimeric aggregates, genotypes reaping group benefits -i.e. 'cheaters'-are expected to bear a selective advantage relative to more 'cooperative' types. In the absence of positive assortment, selection favours cheaters over cooperators, as reproductive gains are amplified across cycles of multicellular aggregation. Such 'tragedy of the commons' scenario [START_REF] Hardin | The Tragedy of the Commons[END_REF][START_REF] Rankin | The tragedy of the commons in evolutionary biology[END_REF] is commonly invoked as a threat to the evolutionary stability of multicellularity, and of particular relevance close to the transition from single cells to higher levels of integration [START_REF] Clarke | Origins of evolutionary transitions[END_REF][START_REF] Rainey | Resolving Conflicts During the Evolutionary Transition to Multicellular Life[END_REF].

In the absence of mechanisms that efficiently purge invading cheaters -such as single-cell bottlenecks in clonally growing organisms -aggregative life cycles are particularly fragile to cheating [START_REF] Strassmann | Altruism and social cheating in the social amoeba Dictyostelium discoideum[END_REF][START_REF] Vos | Social Conflict in Centimeter-and Global-Scale Populations of the Bacterium Myxococcus xanthus[END_REF]. Theory indeed predicts that a sufficient degree of positive assortment, or high relatedness, is necessary for cooperation to be evolutionarily stable [START_REF] Fletcher | A simple and general explanation for the evolution of altruism[END_REF][START_REF] Queller | Genetic relatedness in viscous populations[END_REF]. If cooperators recognize cooperators, for instance through mechanisms of kin recognition or green beard effects [START_REF] Hamilton | The genetical evolution of social behaviour[END_REF], then cooperative groups can be advantaged over excluded cheaters. Even though recognition mechanisms have been identified in Dictyostelium, high segregation efficiency is only attained among distant species [START_REF] Sathe | Interspecies and intraspecies interactions in social amoebae[END_REF]. Moreover, it is unlikely that highly sophisticated recognition systems were present before the multicellular organization was stabilized, leaving the question open of how positive assortment was first established in primitive multicellular life cycles [START_REF] Clarke | Origins of evolutionary transitions[END_REF]. More relevant to understanding the independent repeated evolution of aggregative life cycles [START_REF] Grosberg | The Evolution of Multicellularity: A Minor Major Transition?[END_REF] seem thus to be passive mechanisms of assortment. Ancestral heterogeneity in individual cell properties, indeed, can be the source of positive assortment also in the absence of cell-cell signalling and recognition.

Theoretical explorations have for instance shown that costly adhesive traits are selected if they enhance group cohesion, even when cells reshuffling at every generation brakes identity of descent in multicellular groups [START_REF] Garcia | Group formation and the evolution of sociality[END_REF][START_REF] Garcia | Differential adhesion between moving particles as a mechanism for the evolution of social groups[END_REF][START_REF] Garcia | The evolution of adhesiveness as a social adaptation[END_REF][START_REF] Van Gestel | Phenotypic Heterogeneity and the Evolution of Bacterial Life Cycles[END_REF].

Here, we focus on heterogeneity in motility as a possible source of assortment when multicellular aggregates are generated by self-organization of moving cells, represented as self-propelled disks that can stick to one another upon encounter [START_REF] Szabó | Phase transition in the collective migration of tissue cells: Experiment and model[END_REF]. Heterogeneity in cell motility is observed both in Myxococcus Xantus [START_REF] Vos | Natural variation of gliding motility in a centimetrescale population of Myxococcus xanthus: Gliding motility in Myxococcus xanthus isolates[END_REF] and D. discoideum, where it has been associated to diverse developmental fates [START_REF] Azhar | Cell cycle phase, cellular Ca2+ and development in Dictyostelium discoideum[END_REF] and to efficient foraging [START_REF] Rossine | Structured foraging of soil predators unveils functional responses to bacterial defenses[END_REF]. Differences in motility, moreover, have been proposed to drive the emergence of life cycles with temporally compartmentalized phases of aggregation and disaggregation [START_REF] Miele | Aggregative cycles evolve as a solution to conflicts in social investment[END_REF].

A feature of such cycles, as well as of aggregative life cycles, is the collective timescale that cells have at their disposal to aggregate. After a specific time span, multicellular aggregates proceed to further developmental steps, or are dispersed.

When different strains co-aggregate, the existence of a succession of wellregulated developmental steps allows to measure their reproductive fitness in terms of the fraction of spores (or cysts) they have produced at the moment of dispersal. In Dictyostelium, 'spore bias' of a given strain, i.e. the deviation in the proportion of spores of that strain from its proportion in the initial mix, is thus commonly used as a metric for pair-wise comparison of the social behaviour of different cellular populations [START_REF] Strassmann | Altruism and social cheating in the social amoeba Dictyostelium discoideum[END_REF]Fortunato et al., 2003a;[START_REF] Buttery | Complex genotype interactions influence social fitness during the developmental phase of the social amoeba Dictyostelium discoideum[END_REF][START_REF] Kuzdzal-Fick | An invitation to die: initiators of sociality in a social amoeba become selfish spores[END_REF], 2011;[START_REF] Gilbert | High relatedness maintains multicellular cooperation in a social amoeba by controlling cheater mutants[END_REF]. If aggregates all had a similar composition (no segregation), and cells of different types behaved identically during development, spore bias would be essentially set by the proportion of cells that did not join the aggregates, the so-called 'loners' in D. discoideum [START_REF] Dubravcic | An evolutionarily significant unicellular strategy in response to starvation in Dictyostelium social amoebae[END_REF][START_REF] Tarnita | Evolutionary construction by staying together and coming together[END_REF][START_REF] Martínez-García | Lack of Ecological and Life History Context Can Create the Illusion of Social Interactions in Dictyostelium discoideum[END_REF] or 'peripheral rods' in M. xanthus [START_REF] O'connor | Behavior of peripheral rods and their role in the life cycle of myxococcus xanthus[END_REF]. The propensity of aggregative species to give rise to non-aggregated cells has been proposed to be key to survival in unpredictably fluctuating environments, and to the maintenance of diversity in natural settings [START_REF] O'connor | Behavior of peripheral rods and their role in the life cycle of myxococcus xanthus[END_REF][START_REF] Dubravcic | An evolutionarily significant unicellular strategy in response to starvation in Dictyostelium social amoebae[END_REF][START_REF] Tarnita | Fitness tradeoffs between spores and nonaggregating cells can explain the coexistence of diverse genotypes in cellular slime molds[END_REF]. If it is known that the propensity of cells to aggregate depends on context (both on the abiotic and biotic environmental conditions [START_REF] Rossine | Eco-evolutionary significance of "loners[END_REF][START_REF] O'connor | Behavior of peripheral rods and their role in the life cycle of myxococcus xanthus[END_REF]), the mechanisms determining how many and what cells are more likely to aggregate are still unknown.

We leverage numerical methods developed in active matter physics to explore how heterogeneity in motility affects aggregation biases. Previous studies on Brownian particles with different activity have focused on dense mixtures, where Motility-Induced Phase Separation gives rise to the coexistence of a gas of isolated particles with aggregates, even in the absence of adhesive forces between particles [START_REF] Stenhammar | Activity-Induced Phase Separation and Self-Assembly in Mixtures of Active and Passive Particles[END_REF][START_REF] Rogel Rodriguez | Phase behaviour and dynamical features of a two-dimensional binary mixture of active/passive spherical particles[END_REF][START_REF] Kolb | Active binary mixtures of fast and slow hard spheres[END_REF]. Within aggregates, moreover, motility differences result in sorting, with further implications in multicellular development (? [START_REF] Beatrici | Cell sorting based on motility differences[END_REF].

We characterize finite-time aggregation patterns in low-density binary mixes of adhesive particles that only differ for their intrinsic motility. Given the potential implications of such motility differences in the evolutionary origin of aggregative multicellular life cycles, we propose a criterion to infer aggregation bias in binary mixes from the properties of homogeneous populations. Moreover, we address the dependence of aggregation biases on the frequency of cells of different motility, showing that non-trivial outcomes are possible, though not all that could be expected to occur. A particularly interesting case is the establishment of frequency-dependent selection associated to the evolutionarily stable coexistence of different motility types.

Methods

Source code

Simulations were performed using the Python language. The Pytorch library [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] for tensor computation allowed us to reduce the computational time while working with a large number of particles (N = 10 4 ).

The source code for the simulation can be found online (https://github.com/ MathieuForget/HeteroSpeed_Aggregation.git).

Identification of aggregates

Aggregates are defined as groups of at least 5 aggregated particles. To be considered as aggregated, a particle has to be connected to at least 4 others particles in the final state of the simulation. Two particles are considered connected if their pairwise distance is lower than R 0 .

Variance in aggregates composition

In order to quantify the level of segregation between the two particles types in a binary mix, the variance in aggregate composition was scored and then normalized by its maximum theoretical value. For each simulation, this maximum value is computed by estimating the variance in aggregates composition after having artificially segregated aggregated particles by dividing each type in aggregates of size equal to the mean observed aggregates size. This corresponds to a case where the same number of particles have aggregated, but they only group with those of their same type. The normalized variance in aggregate composition thus goes from 0, when the particles are well mixed in the aggregates, to 1, when the two particles types are completely segregated into distinct aggregates, each composed by only one type of particles (Fig. 3B).

Bias in aggregates composition

At the end of each simulation, the mean proportion of each population in the aggregates is estimated. The bias in aggregate composition of the focal population is measured as the deviation between its mean proportion in the aggregates and its proportion in the mix.

Particles connectivity

The number of particles that are closer than R 0 to a focal particle is used as a proxy for its connectivity. We estimated populations mean connectivity by averaging the connectivity of aggregated particles of a given type. In binary mixes, the relative connectivity of this type was estimated as the deviation between its mean connectivity and that of the other type of particles.

Results

A self-propelled particles model for the aggregation of motile cells.

In order for aggregation to bring together cells of different populations, cells need to be motile. Different aggregative organisms have different types of motility and different shapes. Here, we consider a simple representation of moving cells as persistent Brownian walkers. This representation corresponds to observations of single-cell motility in Dictyostelium discoideum [START_REF] Golé | A quorum-sensing factor in vegetative Dictyostelium discoideum cells revealed by quantitative migration analysis[END_REF][START_REF] Li | Dicty dynamics': Dictyostelium motility as persistent random motion[END_REF], and was used to model collective behaviour of keratocytes populations [START_REF] Szabó | Phase transition in the collective migration of tissue cells: Experiment and model[END_REF]. Single cells are represented as self-propelled disks of equilibrium radius R eq /2 (R eq being the equilibrium distance between two interacting particles) that interact within a finite range R 0 . Interaction forces determine, together with stochastic fluctuations, the position of each disk and the orientation of each particle's velocity. In order to avoid to define their behaviour at the boundaries, particles are let move in a square of side L with periodic boundary conditions.

The equation of motion of a particle i with position r i ∈ [L, L] is:

dr i dt = v 0 n i + µ N j=1 F(d ij ). ( 1
)
where v 0 is the self-propelled speed of particle i, whose direction is given by a polarization unitary vector n i , N is the population size, and F is the interaction force between particles i and j, that only depends on their Euclidean distance

d ij = ∥r j -r i ∥.
The mobility parameter µ ensures dimensional consistence, and is set to 1 without loss of generality.

Cell-cell adhesion acts only on cells that are not too far apart. The model considers a radial, finite-range interaction force F that is piecewise-linear in the distance d ij (Fig. 1A). If d ij ≤ R eq , the two particles i and j repulse each other strongly, in order to account for cell volume exclusion. If R eq < d ij ≤ R 0 , the two particles attract each other with an elastic force of intensity F adh that quantifies adhesion strength. This choice reflects the possibility that cell shape deformation allows cells to be more distant than twice their equilibrium radius.

The radial force is thus defined as:

F(d ij ) = e ij          F rep dij -Req Req d ij < R eq F adh dij -Req R0-Req R eq ≤ d ij ≤ R 0 0 R 0 ≤ d ij (2)
where e ij = ri-rj ∥ri-rj ∥ is the unitary vector connecting the centers of the two particles.

The force F and its action on particles' self-velocity are represented in Figs.

1A and B. In the simulations, time is discretized in steps of (short) duration dt.

Between two successive time steps, interaction forces induce a small incremental change in the angle θ n i (t) = arcsin(n i (t)) of particle self-propelling velocity.

Even though there is no direct velocity alignment in this model, the action of the force is to reduce the deviation of the particle velocity from the scattering direction v i = dri dt defined by eq. 2. Such a relaxation happens on the timescale τ . The equation of motion for the angle is therefore:

dθ n i dt = 1 τ arcsin n i × v i |v i | • e z + η ξ i (3)
where e z is a unit vector orthogonal to the plane of motion. The cross product is null when the velocity and the force are aligned, case in which the direction of motion will not be altered by particles interaction.

In isolation (i.e., in the absence of interaction forces F = 0), particles move as persistent random walkers. Typical trajectories, shown in Fig. 1C, are consistent with previous [START_REF] Golé | A quorum-sensing factor in vegetative Dictyostelium discoideum cells revealed by quantitative migration analysis[END_REF][START_REF] Li | Dicty dynamics': Dictyostelium motility as persistent random motion[END_REF] and our own (Fig. 1D) experimental observation of individual Dictyostelium discoideum cells trajectories.

The orientation θ n i of the self-propelled velocity v i = v 0 n i (Eq. 1) is described by a Wiener process [START_REF] Huang | Stochastic and Multiple Wiener Integrals for Gaussian Processes[END_REF], where a white angular noise term ηξ i is added between successive time steps. ξ i (t) is a random variable chosen with a uniform probability on the interval [-2π, 2π], uncorrelated across successive time steps. The noise intensity η defines the characteristic timescale τ r = 1/η 2 of persistence of the random motion. The larger the angular noise, the less persistent is the motion of the particle.

In the following, we will consider situations where τ r is small (large η) compared to the characteristic relaxation time τ , so that aggregates do not display collective motility (in such regime, extensively studied in active matter [START_REF] Grégoire | Moving and staying together without a leader[END_REF][START_REF] Casiulis | Velocity and Speed Correlations in Hamiltonian Flocks[END_REF], groups are less well-defined than in cases where particles form long-lived, quasi-stationary aggregates). [START_REF] Dubravcic | An evolutionarily significant unicellular strategy in response to starvation in Dictyostelium social amoebae[END_REF][START_REF] O'connor | Behavior of peripheral rods and their role in the life cycle of myxococcus xanthus[END_REF].

Finite-time states of homogeneous populations.

We first explored how the system's parameters affect the patterns of aggregation of particles with identical velocities. Particular attention was given to two parameters that can be controlled or at least measured in experiments, and that are known to influence encounter rates among particles: the packing fraction ρ and the particles' self-propulsion speed v 0 .

The packing fraction is the ratio between the areas covered by the particles relative to the total area of the simulation space. It is thus proportional to the population density, which is straightforward to tune in experimental settings, and is known to critically affect collective behaviour in aggregative species [START_REF] Hashimoto | Cell density dependence of the aggregation characteristics of the cellular slime mould dictyostelium discoideum[END_REF]. As motivated in the introduction, the speed of particle displacement v 0 is interesting because heterogeneity in cell motility is known to exist in aggregative species [START_REF] Goury-Sistla | Bimodal distribution of motility and cell fate in Dictyostelium discoideum[END_REF][START_REF] Vos | Natural variation of gliding motility in a centimetrescale population of Myxococcus xanthus: Gliding motility in Myxococcus xanthus isolates[END_REF], and it can be altered by changing growth conditions in the lab [START_REF] Varnum | The developmental regulation of single-cell motility in Dictyostelium discoideum[END_REF]. Particle encounters become more likely when both packing fraction and velocity are increased. However, increasing one or the other is not equivalent.

The self-organized patterns reached after a finite time by a population of identical particles were examined for different values of speed v 0 and packing fraction ρ (Fig. 2). Simulations were run for a fixed duration t f . This choice reflects the fact that cells do not have an infinite time at their disposal to complete aggregation. In biological populations, such a timescale is set by the duration of different phases of the life cycle or, for ancestral unicellular types, by environmental variation. In the final state -analysed whether the system had converged to its asymptotic configuration or not -aggregates were defined as groups of at least 5 particles as described in the Methods. The setting of a size threshold reflects the observation that in Dictyostelium aggregate speed is positively correlated to size [START_REF] Rieu | Direct mechanical force measurements during the migration of dictyostelium slugs using flexible substrata[END_REF][START_REF] Bonner | Why does slug length correlate with speed during migration inDictyostelium discoideum[END_REF], so that small slugs have impaired collective motility, and possibly development. Therefore, their contribution to successive generations is akin to that of solitary cells.

More generally, size-dependence is a shared feature of many collective functions, which require a quorum of composing cells for selection to start acting on them [START_REF] Cornforth | Synergy and group size in microbial cooperation[END_REF].

At very low packing fractions, the system does not display any sizeable aggregate. Indeed, particles do not have the time to overcome the average distance that initially separates them. The same is true for denser populations when either particles are very slow, or they are very fast. In both cases, the final state is gas-like. For intermediate velocities and larger densities, particles form aggregates of different geometry. Figures 2 and S1 represent,overlaid to the heatmap of the fraction of aggregated particles, the final snapshot of simulations and the variation in time of different population-level statistics.

Different possible qualitative outcomes of the aggregation process illustrated there are discussed below.

At low packing fraction ρ and low velocity v 0 , particles self-organize into a large number of relatively small aggregates (Fig. 2, case A). The isotropy of the clusters reflects the homogeneity of the initial particle distribution. At the end of the simulations, a sizeable fraction of particles are found outside the aggregates. As indicated by the monotonous increase in the fraction of aggregated particles over time (Fig. S1A), this is a finite-time effect and every particle would eventually aggregate for higher t f values (Fig. S2). Under these conditions, the fate of each particle therefore depends on its local environment at the beginning of the simulation. This parallels the hypothesis by Rossine et al that local fluctuations in the concentration of a quorum signal explains why some Dictyostelium cells end up as loners at the end of aggregation [START_REF] Rossine | Eco-evolutionary significance of "loners[END_REF].

D E C A B v 0 ρ % aggregated particles
At high particle speed and low packing fraction (Fig. 2, case B), the number of aggregates quickly reaches a plateau and then slowly decreases, whereas aggregate mean size constantly increases over time (Fig. S1 B). Similarly to an

Ostwald ripening process (Ostwald, 1896), particles evaporate (i.e detach) from small groups to join bigger and more stable aggregates (Fig. S2). Aggregate accretion is essentially limited by particles evaporation, causing a fraction of particles to be always found outside the aggregates.

By continuously increasing particle speed at low packing fraction, the system transitions between the two previously described states where particle clusters coexist with non-aggregated particles. The latter have however different origin depending on particle speed. At low particle speed, aggregation leaves out "latecomer" particles that did not reach any aggregate within the finite time of the simulation; at high speed, the "gas-like phase" is composed by particles that have evaporated from aggregates. Increasing particle speed, the transition between these two scenarios occurs through a phase where aggregation is very efficient. The fraction of aggregated particles, indeed, first increases and then decreases, and reaches a maximum at the optimal speed v opt , which depends weakly on the packing fraction ρ.

At such intermediate velocity and packing fractions (illustrated in Fig. 2, case C), self-organization leads to the formation of multiple rounded aggregates, that contain almost every particle. For much longer simulation times, aggregates would coarsen and progressively merge into one single aggregate of size N , as indicated by the combined decrease in aggregate number and increase in mean size (Fig. S1,C). Such coarsening process however is largely irrelevant for aggregation over finite times.

Finally, at higher packing fraction, particles rapidly self-organize into dendritic structures (Fig. 2, cases D and E). Similarly to what can be observed at lower packing fractions, the few non-aggregated particles at low speed are "latecomers", whereas at high particle speed, non-aggregated particles are "gas-like particles", where -as for MIPS -steric interactions limit the coexistence of multiple independent aggregates [START_REF] Fily | Athermal Phase Separation of Self-Propelled Particles with No Alignment[END_REF]. The measure of aggregates size is not meaningful under these conditions, the packing fraction being so high that every particle is connected to a large number of other particles virtually from the start.

Differences in motility bias the outcome of aggregation Let us first consider two 1:1 binary mixes where particles have self-propelled velocities that are either both to the left, or both to the right of v opt , so that non-aggregated particles stem from the same of the two previously discussed mechanisms (Fig. 3A).

Aggregation biases can manifest as differential fractions of particles that aggregate, in aggregate composition, or in the topology of interactions within aggregates. In order to separately discuss those features, we quantified three statistical properties or order parameters (illustrated in Fig. 3B): segregation, aggregation bias and connectivity. These metrics, moreover, are easy enough to be measured in real populations, as long as there is a way to visually distinguish two co-aggregating populations.

Segregation measures the degree of heterogeneity in aggregate composition.

Under high segregation, which could be expected if particles assorted preferentially with those of the same type, the overall behaviour would be the superposition of the patterns of each composing population in isolation. By analysing the final snapshot of a simulation, every aggregate was characterized by its composition (fraction of particles of each type). The variance of this quantity was normalized by the maximal theoretical variance given the observed number of aggregates and their size, computed when each aggregate is composed of only one type of particles (see Methods). A null value of this metric indicates that the two types are represented in similar proportion in all aggregates, whereas this parameter converges towards 1 as the two types separate in different aggregates (Fig 3B). In biological chimeras where cells are able to discriminate and preferentially aggregate with kins [START_REF] Gruenheit | A polychromatic 'greenbeard' locus determines patterns of cooperation in a social amoeba[END_REF], this order parameter may therefore be close to 1, whereas in non-segregating populations it is expected to have small values (only vanishing in the limit of infinite population sizes, where fluctuations in the initial position of particles become irrelevant).

Aggregation bias was then estimated as the deviation between the mean aggregate composition and the proportion of each type of particles in the initial mix. This metric indicates to what extent differences in particle speed skew their ability to aggregate. In aggregative microbes only cells that manage to partake multicellular development have a chance to survive prolonged starvation. The aggregation bias measured in chimeras can thus be used as a proxy to connect the outcome of one cycle of aggregation with the long term success of a population, in the absence of further developmentally-induced divergence in cell fate [START_REF] Dubravcic | An evolutionarily significant unicellular strategy in response to starvation in Dictyostelium social amoebae[END_REF][START_REF] Tarnita | Evolutionary construction by staying together and coming together[END_REF].

Finally, we quantified particles mean within-aggregate connectivity. For each type of aggregated particles, connectivity is measured as the number of particles (of any type) to which a focal particle is connected within an aggregate.

Two particles are considered connected at a given time if they are closer than the maximal distance of interaction R 0 (see Methods). High connectivity is associated to a central position within an aggregate, whereas particles with a lower connectivity are located closer to aggregates periphery. Connectivity thus reflects different spatial arrangements within clusters. In extant aggregative microbes, a more central or peripheral positioning can be associated to differential exposition to stress [START_REF] Smukalla | FLO1 is a variable green beard gene that drives biofilm-like cooperation in budding yeast[END_REF] or to differentiation cues [START_REF] Julien | Spatial control of cell differentiation in imyxococcus xanthus/i[END_REF]. We observed that the level of segregation is comparably low in the two mixes (Fig. 3C), suggesting that different aggregates have similar composition, bar statistical fluctuations. In other words, heterogeneity in particles motility does not cause segregation of the two types and the mean aggregate composition is representative of every group. It is therefore meaningful to use it to quantify the effect of heterogeneity in particle speed on the outcome of aggregation.

By analysing biases in mean aggregates composition relative to the 1:1 composition of the mix, we found that the slower type was under-represented in the aggregates (i.e. with a negative bias) in the DA mix, whereas the opposite results was found in E mix. Heterogeneity in particle speed thus biases the representation of the two co-aggregating populations in the aggregates, so that particles that aggregate more efficiently in isolation are over-represented in the aggregates. These are the type whose fraction of aggregated particles at t = t f is higher (Fig. 3A), i.e. the faster type in the DA mix and the slower type in the E mix.

Mean particle connectivity of the two particles types indicates that, moreover, particles that aggregate more efficiently in isolation are enriched in the core of the aggregates. The type whose speed is closer to v opt seems to nucleate the aggregation of the other type and ends up being over-represented in the aggregates at the end of the simulation.

Differences in motility result into three types of frequency-dependent aggregation biases.

The process of self-organization responsible for the aggregation biases discussed above results from interactions between particles with different speed. It is hence expected to be sensitive to the proportion of each type of particles in the mix, which underpin the probability that particles of the same or different types encounter. Frequency-dependent biases have been moreover observed in aggregative multicellular species: the intensity of spore bias in Dictyostelium, and sometimes even its sign [START_REF] Sathe | Complex interactions underpin social behaviour in Dictyostelium giganteum[END_REF][START_REF] Madgwick | Strategic investment explains patterns of cooperation and cheating in a microbe[END_REF], depend on the proportion of cells in the binary mix. Such frequencydependent social interactions can lead, on the evolutionary timescale, to regimes where contrasting social strategies coexist, as we will discuss later. In order to identify the possible paths open to the evolution of aggregation in populations with heterogeneous cell motility, we investigated the effect of binary mix composition on aggregation bias. Let f 1 be the proportion of particles with speed v 1 in the binary mix (then

f 2 = 1 -f 1
is the proportion of particles with speed v 2 ). As indicated in Fig. 4, there are only four qualitative types of biases that are continuous in the frequency f 1 and change sign at most once in the interval [0, 1]. These scenarios can be classified by looking at the variation of the bias for extreme values of f 1 : when the slope of the bias has the same sign at low and high frequency, the type associated with negative bias gets excluded over evolutionary times (Fig. 4A andB); when there is positive (negative) frequency-dependence at low (high) frequency, coexistence of the two types is evolutionary stable (Fig. 4C); when there is negative (positive) frequency-dependence at low (high) frequency, the monomorphic populations are two alternative evolutionary stable solutions (Fig. 4D). Assuming that every aggregated particle has the same probability of being transmitted to the next generation, which is reasonable for primitive aggregative life-cycles, the shape of the bias in aggregate composition determines the outcome of evolution over repeated cycles of aggregation and dispersal, as exemplified in Fig. 4A.

We examined what types of bias could be observed in binary mixtures of particles of different speed, and tried to establish a link to the properties of homogeneous populations of each type on its own. Aggregate composition was analysed in binary mixes of different composition, keeping the packing fraction equal to 0.07. In order to reduce simulation time, we considered only extreme f 1 values (0 < f 1 < 0.05 and 0.95 < f 1 < 1).

We expect the aggregation process to be influenced both by the absolute and by the relative value of the velocities of the two types of particles. We have thus explored separately two scenarios: first, the speed of the first type was fixed and that of the second type was changed; second, both speeds were changed, but their difference was kept constant.

In the first analysis, the speed of particles of the focal population was set to the 'optimal' value (v 1 = v opt = 13) and the speed differential ∆v = v 1v 2 between the two particles types was gradually changed from negative to positive. This scan (Fig. S3) confirms the intuition that the better aggregator is always over-represented in the aggregates when mixed with particles with a different speed, whatever the mix composition f 1 . Indeed, the bias in aggregate composition of the better-aggregator (i.e. the v 1 -population) is systematically positive at both low and high f 1 values. Moreover, the slope of the bias for low and high f 1 increases with |∆v|, indicating that the bias intensity is positively 20 correlated with the level of heterogeneity between the two particles types.

In the second type of binary mix, where the speed differential ∆v is kept constant while particle speeds are simultaneously increased, the slopes at extreme frequencies instead change sign (Fig. S4). Besides the two cases where one particle type dominates for any frequency (Fig. 4 A and B), we observed also the frequency-dependent profile that predicts long-term coexistence (Fig. 4 C). The latter occurs when binary mixes are composed of particles whose speed is on opposite sides of v opt (Fig. S4, third and fourth row). The bias is in favor of the slower particles when they are at low frequency in the mix, whereas at high frequencies they are under-represented in the aggregates. The shaded area represents the 95% confidence interval.

Biases in aggregates composition however appear rather weak and noisy (as indicated by the large confidence intervals). We thus ran 20 independent simulations for f 1 values ranging from 0 to 1 and using a larger speed differential (∆v = 6) to confirm that binary mixes composed of particles whose speed falls on different sides of v opt sustain long-term coexistence. As expected, this combination yielded a more conspicuous and less variable bias (Fig. 5B).

Mean connectivity of aggregated particles for the two particles types provided mechanistic information on the underpinning of such frequency-dependent bias. Figure 5D shows that particles with higher connectivity -i.e. more central in the aggregates -are over-represented in the aggregated fraction, irrespective of their speed, consistent with previous observations (Fig. 3C). When slow particles are in low proportion in the mix, aggregates are essentially composed of fast particles. In such aggregates, fast particles have lower connectivity (Fig. 5C and D), i.e a peripheral location, so that they can be pushed out of the aggregate by stochastic fluctuations and evaporate. Slow particles, on the other hand, have a high chance of being integrated in an aggregate quickly formed by fast particles, hence the positive bias in aggregates composition is in their favor.

When slow particles are in high proportion in the binary mix, they prevail in the aggregates. Due to their higher speed, fast particles join aggregates earlier and are thus located closer to the aggregates center than slow particles (Fig. 5C).

The central location of fast particles prevents their evaporation, and delayed aggregation simultaneously causes slow particles to be under-represented in the aggregates. On longer timescales, however, all slow particles join an aggregate, and fast particles are expected to sort to the periphery of aggregates (?). The bias against particles that are too slow to aggregate is therefore expected to progressively vanish as the time t f available to aggregate increases.

Numerical simulations of the binary mix of Fig. 5 were run for twice and five times longer time spans (t f = 200 and t f = 500). As predicted, the range of f 1 values for which the slower particles are under-represented in the aggregated fraction decreases accordingly, and vanishes for the longest period (Fig. 6).

When aggregation lasts long enough, so that all particles end up aggregating, evaporation of fast particles is thus the only source of aggregation biases. Hence, only the bias in favour of slower particles (Fig. 4) is asymptotically observed.

As a consequence, the existence of life cycles that bound the duration of the aggregation phase could support the evolutionary maintenance of highly motile cells, being otherwise counter-selected, and simultaneously the establishment of a polymorphic equilibrium.

Alongside particle speed, adhesion strength is certainly a parameter that plays a key role in determining the outcome of an aggregation process. Moreover, it is likely to be heterogeneous in cellular populations. The complete examination of the interplay of speed and adhesion heterogeneity is beyond the scope of this work, however it is important to consider the extent to which the conclusions drawn so far keep holding when adhesion strength is varied.

Biases induced by heterogeneity in particles motility depend predictably on particles adhesion strength.

In addition to particle speed and packing fraction (Fig. 2), particle adhesion strength is also expected to affect aggregation efficiency. This parameter indeed determines the probability that two particles stick together when they encounter. Increasing particles adhesion strength can thus be expected to reduce the rate at which particles evaporate from aggregates, hence affect the balance between evaporation and delayed aggregation.

The fraction of aggregated particles at the end of the simulation was estimated for a range of particle speeds and adhesion strengths in a population with homogeneous speed in order to obtain a quantitative understanding of how particle adhesive strength influences the outcome of aggregation. Figure S5 displays the proportion of aggregated particles for packing fraction ρ = 0.07.

Variation in F adh does not affect the outcome of aggregation for small particle speeds. In this regime, aggregation over a finite time is indeed limited by the low encounter rate of isolated slow particles. As soon as particles are fast enough to interact, almost any adhesion strength except the weakest leads to aggregation.

The fraction of aggregated particles increases with their speed to the point where evaporation starts being important. The higher v 0 , the higher the force intensity required to counterbalance the effect of evaporation. In the Supplementary Information, we present a microscopic model that predicts the quadratic relation between F adh and particle speed observed in the simulations (Fig. 7A).

Let us now consider a binary mix of particles with different speed, where adhesion strength varies uniformly for both particle types, while other parameters are the same as for Fig. 5 (v 1 = 10, v 2 = 16, ρ = 0.07). We previously proposed that the type of bias in aggregates composition depends on the relative position of v 1 and v 2 , relative to v opt . We can now check that this assumption correctly predicts the change of the bias when v opt changes as a function of F adh .

For any fixed couple of velocities v 1 and v 2 , there are three possible ranges of F adh where v opt is located to the left, to the right, or in between their values, as illustrated in Fig. 7. By choosing adhesion strengths that belong to such Overall, when the two co-aggregating populations are on both side of v opt , aggregates composition is biased because of the interplay of two counteracting processes: evaporation of fast particles and delayed aggregation of slow parti-cles. The relative contribution of these two factors determines the sign of the bias and depends on the balance of the types in the binary mix. When the speed of both particle types is smaller than v opt , the slower is under-represented in the aggregates due to delayed aggregation. When both speeds are larger than v opt , the slower type is over-represented in the aggregates due to the higher evaporation of faster particles. Finally, when v 1 < v opt < v 2 , the competition between delayed aggregation of the slower particles and evaporation of the faster particles provides different results if one or the other particle is at low frequency, giving rise to a frequency-dependent bias. As long as the duration of the aggregation phase is finite, then, the evolutionary dynamics is predicted to stabilize the coexistence of different types of particles. Interestingly, the frequency-dependent bias associated to bistability (Fig. 4D) was not observed. This scenario would require particles faster than v opt to be over-represented in the aggregates when abundant (low f 1 values), which can only be achieved if differential evaporation in favor of fast-particles can be compensated. Theoretically, this may occur if fast particles were surrounded by a "shell" of slow particles, which, for initially uniformly assorted particles, would require additional mechanisms of spatial sorting. It may, for instance, be achievable if slower particles were also less adhesive [START_REF] Garcia | Differential adhesion between moving particles as a mechanism for the evolution of social groups[END_REF].

These results show that, somehow unexpectedly given the complexity of the self-organized dynamics and the transient nature of aggregation patterns, biases in binary mixes can be predicted once the behaviour of the two populations in isolation is characterized.

Discussion

When multicellularity emerged from unicellular ancestors, a transition occurred from populations under selection for individual properties to a higher level of biological organization, where collective adaptation became possible.

Aggregative multicellularity, one of the two routes followed in this transition, involves the integration of cell aggregation into a life-cycle that also encompasses free living. Cell aggregation offers a mechanism of fast and transient grouping, that may be adaptive for coping with rapidly and unpredictably changing environments [START_REF] Márquez-Zacarías | Why have aggregative multicellular organisms stayed simple?[END_REF]. As a side-effect of fast grouping, multicellular aggregates are however generally composed of a mosaic of cells in different phenotypic states. Also observed in extant aggregative species (Fortunato et al., 2003b;[START_REF] Sathe | Genetic Heterogeneity in Wild Isolates of Cellular Slime Mold Social Groups[END_REF][START_REF] Gilbert | High relatedness maintains multicellular cooperation in a social amoeba by controlling cheater mutants[END_REF], single-cell heterogeneity between co-aggregating cells may have been particularly important for the early steps of the transition to aggregative multicellularity, as it can promote the recruitment of ancestral functions towards incipient division of labour. On the other hand, before mechanisms of assortment evolved to reduce it, intra-group heterogeneity might have caused widespread genetic conflicts that threatened multicellular organization.

Differences in cell properties between co-aggregating cells is known to bias their probability to differentiate into dormant phenotypes at the end of multicellular development [START_REF] Forget | Social conflicts in Dictyostelium discoideum : a matter of scales[END_REF]. However, understanding the evolutionary implications of such single-cell heterogeneity on the outcome of aggregation is currently limited by a lack of mechanistic explanations, in particular of how initial diversity translates into biases in the composition of multicellular aggregates. Among the many features that display cell-to-cell variation, cell motility can be expected to directly impact the probability that a cell joins a multicellular aggregate, and therefore affect spore bias. Indeed, as already pointed out for D. discoideum loners [START_REF] Dubravcic | An evolutionarily significant unicellular strategy in response to starvation in Dictyostelium social amoebae[END_REF][START_REF] Tarnita | Fitness tradeoffs between spores and nonaggregating cells can explain the coexistence of diverse genotypes in cellular slime molds[END_REF] only aggregated cells can differentiate into spores and survive the dispersal phase of the life cycle.

We used individual-based simulations to confirm the intuition that heterogeneity in cell motility affects the outcome of aggregation. Heterogeneity in single-cell properties is known to act within multicellular aggregates, where cell sorting can be a source of bias in the developmental fate of cells [START_REF] Jones | Relation between the rate of cell movement under agarose and the positioning of cells in heterotypic aggregates[END_REF][START_REF] Houle | A glycosylation mutation affects cell fate in chimeras of Dictyostelium discoideum[END_REF][START_REF] Umeda | Cell sorting by differential cell motility: a model for pattern formation in Dictyostelium[END_REF][START_REF] Strandkvist | A kinetic mechanism for cell sorting based on local variations in cell motility[END_REF]. We showed that biases can also emerge during aggregation, before cells experience proper 'social' interactions. In primitive aggregative life-cycles, before the evo-lution of cell coordination mechanisms susceptible of reshuffling or canceling out aggregation biases, cell motility properties inherited from a unicellular ancestor may therefore have been recruited for their effects on group adaptation.

Numerical exploration of different binary mixes highlighted three types of aggregation biases resulting from heterogeneity in particles motility. Notably, when particles that are too slow to join a group before the end of aggregation are mixed with particles that are fast enough to leave an aggregate, the sign of the bias in the aggregated fraction changes with frequency. Slow particles are over-represented in the aggregates when they are rare in the population, but the opposite is true when they are in a high proportion. This type of frequencydependent bias allows to stabilize polymorphism in populations of interacting genotypes [START_REF] Hudson | Altruism, Cheating, and Anticheater Adaptations in Cellular Slime Molds[END_REF][START_REF] Brännström | Evolutionary dynamics of altruism and cheating among social amoebas[END_REF]. Frequencydependent spore bias has been invoked [START_REF] Madgwick | Strategic investment explains patterns of cooperation and cheating in a microbe[END_REF][START_REF] Parkinson | A Simple Mechanism for Complex Social Behavior[END_REF] as a possible mechanisms for the maintenance of the diversity of strains observed in natural populations of aggregative species, such as Dictyostelium [START_REF] Francis | Genetic structure of a natural population of Dictyostelium discoideum, a cellular slime mould[END_REF]Fortunato et al., 2003b) or Myxococcus xanthus [START_REF] Kraemer | Endemic social diversity within natural kin groups of a cooperative bacterium[END_REF][START_REF] Vos | Social Conflict in Centimeter-and Global-Scale Populations of the Bacterium Myxococcus xanthus[END_REF]. However, it is typically associated to developmental biases that emerge in the course of the multicellular stage of the life cycle( [START_REF] Madgwick | Strategic investment explains patterns of cooperation and cheating in a microbe[END_REF][START_REF] Parkinson | A Simple Mechanism for Complex Social Behavior[END_REF].We have shown that such biases can also emerge passively from heterogeneity in singlecell properties during aggregation. As such, they may have stabilized genetic chimerism early in the transition to aggregative multicellularity, even before development became highly regulated and contributed to limit genetic conflicts in primitive aggregative life-cycle.

Beyond its constructive role on fostering polymorphism during aggregation, heterogeneity in cell speed can also play an active role in later developmental stages. Motility-induced cell sorting within multicellular aggregates [START_REF] Jones | Relation between the rate of cell movement under agarose and the positioning of cells in heterotypic aggregates[END_REF][START_REF] Beatrici | Cell sorting based on motility differences[END_REF] indeed paves the way towards the exploitation of morphogen gradients for establishing differential cell fates. Hence, biases in aggregates composition induced by heterogeneity in cell speed during aggregation could have concomitantly fueled the evolution of division of labor. In extant aggregative species such as D. discoideum and M. xanthus, position within an aggregate is known to determine cell probability to form a spore [START_REF] Araki | Symmetry breaking in Dictyostelium morphogenesis: evidence that a combination of cell cycle stage and positional information dictates cell fate[END_REF][START_REF] Huang | Cells at the Center ofDictyosteliumAggregates Become Spores[END_REF][START_REF] Julien | Spatial control of cell differentiation in imyxococcus xanthus/i[END_REF]. In a model describing cell arrangement within Dictyostelium slugs, slower cells sort to the back during migration [START_REF] Umeda | Cell sorting by differential cell motility: a model for pattern formation in Dictyostelium[END_REF], thus populating preferentially the region that will eventually turn into the spore mass. Developmental sorting can act in the same or opposite direction relative to aggregation biases.

In the evaporation-limited region, faster particles are doubly disadvantaged, as they would be both over-represented in the loners and in the stalk. Motility differences instead act in opposite directions, and could even cancel out, when non-aggregated particles are mostly slow particles delaying their aggregation. In this region, slow particles are indeed under-represented in the aggregates, but once they have joined an aggregate, they are expected to sort to the back of the slug hence to differentiate into spores. Combining aggregation and multicellular development delineates an evolutionary trade-off between being sufficiently fast to join an aggregate during aggregation on the one hand, and sufficiently slow to sort to the back of the slugs on the other hand.

When cell motility is a heritable trait, aggregation biases induced by heterogeneity in single-cell properties underpin its evolutionary trajectory over multiple cycles of aggregation and dispersal. We can imagine that cell motility evolves by the appearance of mutant cells whose speed is a small variation of the resident population. The bias in aggregate formation then determines the likelihood that the mutant trait substitutes the resident, as illustrated in Fig. 4.

Our results indicate that, as long as both the mutant and the resident traits are on the same side of v opt , selection will push particles to approach the optimal, intermediate velocity for which aggregation is maximally efficient. Indeed, the better aggregator in isolation is always over-represented in the aggregates, so that evolution would climb the gradient of the proportion of aggregated particles, that is concave and has a unique maximum (Fig. 3A). However, if the resident and the mutant fall on the opposite sides of v opt -which is increasingly likely, as the optimal value is approached -polymorphic solutions become avail-able, where particles with different motility can coexist. Evolutionary branching is possible in such a situation, as predicted in models for the adaptive dynamics of social strategies [START_REF] Doebeli | The Evolutionary Origin of Cooperators and Defectors[END_REF]. As discussed above, and argued elsewhere [START_REF] Miele | Aggregative cycles evolve as a solution to conflicts in social investment[END_REF], coexistence of slow and fast cells may be integral to the establishment of aggregative life cycles, perhaps even before the two phenotypes were unified as plastic manifestations of one single genotype.

During the last decades, emergent patterns of populations of self-propelled particles have been at the core of research in active matter physics [START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF]. In particular, motility-induced phase separation has been evidenced as a state where steric interactions generate, in sufficiently dense populations, the coexistence of gas and aggregated phases [START_REF] Fily | Athermal Phase Separation of Self-Propelled Particles with No Alignment[END_REF].

We focused instead on smaller densities, where we could observe the emergence of multiple particle aggregates, similar to a field of aggregation of slime moulds.

At these packing fractions, adhesion plays a fundamental role in determining the aggregation pattern. The probability that cells remain together upon encounter indeed derives by a balance between adhesive forces and self-velocity, modulated by fluctuations. Several qualitatively different particle arrangements have been identified as steady-states of populations whose particles are both adhesive and align their motion [START_REF] Grégoire | Moving and staying together without a leader[END_REF]. Before stabilizing their distribution on asymptotic patterns, active particles also display transient clustering (Ostwald, 1896;[START_REF] Van Der Linden | Interrupted motility induced phase separation in aligning active colloids[END_REF]. Such long-lasting regimes can be relevant for biological populations, where aggregation occurs over a finite timescale. Including non-stationary patterns led us to identify two possible regimes where clusters coexist with isolated particles. At low speed, non-aggregated particles are "latecomers" that did not have the chance to interact with other particles by the end of the simulation. Their proportion in the aggregated fraction depends on the duration of aggregation, other than on particle density, and is expected to converge to one as the system reaches its stationary state. A consequence of time-constrained aggregation is that, when they are rare, slow particles are overrepresented in aggregates relative to particles fast enough to evaporate, but have a negative bias when common. Such a frequency-dependent bias is not found in simple models where diversity among interacting agents is depicted in terms of fixed social strategies. Game theoretical representations of interactions within multicellular bodies indeed typically predict that 'cheater' types, that have by definition a positive bias, would be evolutionary stable in the absence of mechanisms favouring interactions among kins [START_REF] Strassmann | Altruism and social cheating in the social amoeba Dictyostelium discoideum[END_REF][START_REF] Strassmann | Evolution of cooperation and control of cheating in a social microbe[END_REF]. Contrary to explanations based on strategic decision-making of cells during multicellular development [START_REF] Madgwick | Strategic investment explains patterns of cooperation and cheating in a microbe[END_REF], we propose here that frequency-dependent biases can also be established by passive mechanisms, based on individual heterogeneity. It will be the subject of future studies whether such dependence is robust to heterogeneity in other parameters. In particular, adhesion on top of velocity could be chosen as a heterogeneous trait, possibly correlated with speed.

These extensions can be informed by experimental observations of aggregative multicellular microbes, where parameters -at least those controlling singlecell motility may be estimated for specific populations. The possibility that frequency-dependent biases are induced during the aggregation phase could be tested in extant aggregative microbes.The aggregation dynamics of populations with different motility could be assessed when they aggregate separately, and in chimeras with different proportions of two types. We are not aware of noninvasive methods to modulate Dictyostelium motility. One possible strategy to generate populations with differential motility is to apply contrasting selective pressures in a directed evolution experiment [START_REF] Adiba | Evolving social behavior through selection of single-cell adhesion in dictyostelium discoideum[END_REF].

More generally, a better understanding of the mechanical properties of unicellular organisms may prove essential to capture the diversity of paths leading to current forms of multicellular organization. to the rapid formation of a large number of small aggregates. On longer timescales, however, particles evaporate from small and unstable aggregates and attach to larger-size aggregates that are more stable (A, snapshots of the system state at t = 50, t = 150 and t = 250).

Conclusion

Similarly to Ostwald's ripening processes (Ostwald, 1896), this results in a decrease in the number of aggregates (B). The fraction of aggregated particles (C) saturates at a value smaller than 1, indicating that in the asymptotic state the aggregated phase coexists with the gas phase, as in MIPS [START_REF] Fily | Athermal Phase Separation of Self-Propelled Particles with No Alignment[END_REF]. In this simulation, ρ = 0.07, v 0 = 17.5, and other parameters were set as indicated in table S1. R 0 /v 0 = 0.01. This value is sufficiently small so that, between two successive time-steps, two particles cannot cross each other's interaction range without interacting.

The interaction parameters F rep and F adh were set so that particles do not collapse (volume exclusion) and can detach from aggregates (evaporation) when v 0 is increased. The equilibrium radius R eq and the interaction range R 0 are such that particles only interact with their direct neighbours, in order to model cell interaction through physical contact. Interactions mediated by diffusing In order to compute the final time of numerical integration t f that corresponds to the finite duration of aggregation, we used the following relationship obtained by dimensional analysis:

t f = t exp f v exp R exp eq
R eq v 0 ≈ 100 (4) Link between optimal velocity v opt and adhesion force F adh Simulations of the self-propelled particles system suggest that v opt and F adh scale like:

v opt = a F adh ( 5
)
where a is a constant (Fig. S5, white dotted line, and Fig. 7A, black line).

In order to understand such a relationship from a theoretical point of view , let us consider the collision between a self-propelled and a stationary particle (for instance, one that belongs to an aggregate -whose average speed is very small), as illustrated in Fig. S6. The equations of motion for each particle in the simulation are:

dr i (t) dt = v 0 n + µF ≡ v (6) dθ(t) dt = 1 τ arcsin(n × v v ) + ηξ j (t) (7)
where r i is the position of particle i, v 0 its self-propelled speed, n a polarization unitary vector with angle θ, µ the mobility, τ is the characteristic time for the polarization to align in the scattering direction defined by v, η is the noise intensity which defines a characteristic timescale τ r ∼ 1/η 2 of the white rotational noise ξ(t). For the analysis that follows, we consider a large value for the alignment and rotational characteristic times, so that the polarization direction (not the scattering direction) remains constant.

We want to estimate the duration of the contact between two particles when they encounter, that is the time they spend at a distance of less than the interaction range R 0 . Comparison between this timescale and the timescale necessary for a particle to escape from an interaction will provide a scaling relation between the force intensity and the optimal velocity.

Consider a particle that is travelling with velocity v approaching the stationary particle. When the incident particle falls within the range R 0 , one has to take into account the deviation in its trajectory that results from the attractive force. Since the force acts radially, we have to evaluate the velocity in the radial direction (defined pointing outwards relative to the stationary particle) resulting from the self-propulsion and the action of the force. This is the projection in that direction of eq. 6:

v r = v 0 cos α + µF adh , ( 8
)
where α is the angle between the radial direction and the direction of the particle velocity (Fig. S6).

In the absence of angular noise, the incident particle will remain within the range of the force as long as the velocity in the outward radial direction is small enough to compensate the last term in eq. 8. The two particles detach when v r is zero, which happens at a critical angle α * where the components of the force and self-propulsion have opposite signs. That occurs when cos α * = µF adh v0 (Fig. S6). A small value for α * implies a small value for the outward velocity component and, consequently, a longer contact time. We may then estimate the contact time as proportional to the attraction force range, R 0 and cos α * and inversely proportional to the speed, v 0 , that is, τ ctc ≈ R0 cos α * v0

. From these relations, one obtains the following scaling:

τ ctc ≈ R 0 µF adh v 2 0 (9)
The larger τ ctc , the longer two particles remain attached.

For a given F adh , let us consider a population of particles with high v 0 so that aggregation is limited by evaporation. Decreasing v 0 results in an increase in τ ctc , particles thus escape at a lower rate and the fraction of aggregated particles increases. When τ ctc > τ , however, whatever the speed, the incident particle will remain bound, as its self-propulsion will be quickly reoriented along a tangent direction. Once the particles enter such an 'orbit' upon contact, further increase in the contact time (corresponding to a decrease in v 0 ) has no effect. A decrease in the fraction of aggregated particles can nonetheless occur because of delayed aggregation, and more so, the slower are the particles. The optimal v 0 (i.e.

v opt ) for which the fraction of aggregated particles is the highest is thus reached when τ ctc = τ .

This condition yields:

v opt = √ R 0 µF adh τ ( 10
)
Figure S6: Timing of particles interaction. A self-propelled particle (blue, represented at successive times) comes into contact with a stationary particle (grey). The size of the stationary particle is enhanced for the sake of visualization. The alignment and rotational characteristic times are considered sufficiently large so that the polarization direction of the self-propelled particle remains constant during the encounter with the stationary one. The self-propelled particle detaches from the stationary particle when the projection of its selfpropulsion in the direction connecting the particle centers is larger than the adhesive strength.

.

Figure 1 :

 1 Figure 1: Self-propelled particle model for migrating and interacting cells. (A) Finite-range interaction forces F(d ij ) between two particles of positions r i and r j and velocities v i and v j (accounting for both auto-propulsion and interaction forces, according to eq. 1) contribute to align n ′ i and n ′ j (n ′ = n(t + dt)), the direction of particles' auto-propelled velocities. These forces only act when the Euclidean distance d ij = |r i -r j | ≤ R 0 . See text for the definition of the other variables. (B) Radial component of the force between the two particles. When d ij < Req, particle i exerts a strong repulsive force on particle j . Otherwise, within the range of interaction R 0 , the latter attracts particle j with a weaker elastic force, representing adhesion between different cells. (C) Typical trajectories of isolated particles during a simulation (t f = 100) using parameters values indicated in table S1 indicate that, in the absence of interactions, particles follow a persistent random walk resembling that of Dictyostelium's single cells. (D) Trajectories of starved Dictyostelium cells recorded every 30 seconds during 1 hour.

Figure 2 :

 2 Figure 2: Finite-time state diagram for a population of identical particles of speed v 0 and packing fraction ρ. The fraction of aggregated particles at the end of the simulation (t f = 100) is illustrated as a heatmap, indicating that the highest fraction of aggregated particles is reached for intermediate particle speeds. Overlaid, snapshots illustrate the final particle distribution for five cases (A, B, C, D, and E), whose aggregation dynamics is discussed in more detail in Supp. Fig. S1. Fixed parameters values are those indicated in table S1. The packing fraction was varied by changing L, keeping N constant.

Figure

  Figure 3C displays these three order parameters for two binary mixes of particles with different speed, as per Fig. 3A. We refer to those as 'delayed aggregation' (DA) mix and 'evaporation' (E) mix to signify the different origin of non-aggregated particles.

Figure 3 :

 3 Figure 3: Effect of heterogeneity in cell motility on the outcome of aggregation in two different types of binary mix. (A) two populations of particles with different motility (DA mix: v 1 =3 and v 2 =6, E mix: v 1 =13 and v 2 =16) where mixed in equal proportion (N = 10 4 , ρ = 0.07 for the total population). (B) Order parameters allowing to capture the effect of heterogeneity in particles motility on the outcome of aggregation. (C) Populations segregation, aggregation bias and mean connectivity within the aggregates estimated at the end of 6 independent simulations as described in B.

Figure 4 :

 4 Figure 4: Four qualitatively different frequency-dependent aggregation biases correspond to distinct evolutionary stable states. Schematic representation of the four simplest classes of bias that can be obtained in a binary mix. The evolutionary dynamics of a binary mix composition (f ) over many rounds of co-aggregation is represented by the yellow arrows. Stable equilibria are represented as filled circles and unstable ones as empty circles. The yellow dotted line in (A) illustrates the construction of the evolutionary trajectory: each round starts with a mix composition equal to the aggregates composition at the end of the previous round, that is obtained by adding to the initial composition the aggregation bias. The four classes of bias give rise to different evolutionary stable solutions: fixation of the focal type (A) exclusion of the focal type (B), coexistence of both types (C) and bistability (D).

Figure 5 :

 5 Figure 5: Heterogeneity in particle speed induces biases whose sign depends on composition of the binary mix. (A) Binary mix of particles with speed on both sides of vopt. (B) Aggregation bias of the slower population when mixed in different proportion (f 1 ) with the faster population. (C) Zoom on snapshots of the final state of the system for f 1 = 0.1 and f 1 = 0.9. (D) Mean connectivity of the two types of aggregated particles as a function of f 1 .(B and D) Twenty independent simulations were run were for each f 1 value.

Figure 6 :

 6 Figure 6: Prolonged aggregation mitigates the frequency-dependent inversion of bias. Aggregation bias of the v 1 -population as a function of its frequency f 1 in a binary mix of particles (v 1 = 10 and v 2 = 16), estimated after different simulation time (t f ). 20 replicates were performed for each f 1 value. The shaded area represents the 95% confidence interval.

Figure 7 :

 7 Figure 7: Particle adhesion affects the outcome of aggregation in binary mixes of particles with heterogeneous speed. (A) Particle speed vopt for which the proportion of aggregated particles at t = t f is highest increases with the intensity of adhesive forces F adh (see Fig. S5). The curve vopt ∼ √ F adh (see Supplementary Information for a derivation from microscopic principles) was fitted to the optimal speed obtained from simulations (grey dots). Three F adh ranges can be defined depending on the relative position of v 1 , v 2 and vopt (indicated by different colours), as previously discussed for varying v 0 . (B) Aggregation bias of the slower population measured at the end of simulations of the binary mix described in Fig. 5, for adhesive strength within the ranges indicated in the left panel (F adh = 4, F adh = 7 and F adh = 12).

Figure S1 :Figure S2 :

 S1S2 Figure S1: Variation in particle speed and packing fraction produces different aggregation dynamics. Snapshots of the final state (first column) of the system and evolution of aggregate number and aggregate mean size (second column), as well as the fraction of aggregated particles over time (third column), for simulations ran using the 5 pairs of parameters values displayed in Fig. 2. The five different regimes discussed in the text correspond to different variations of the three order parameters.

Figure S3 :

 S3 Figure S3: Bias in binary mixes where the focal population has constant speed, and speed differential varies. The first column displays the fraction of aggregated particles when each population aggregates separately at the same total packing fraction ρ = 0.07 (other parameters were set as indicated in table S1). Red dots indicate the constant speed v 1 = 13 of the focal population, blue dots the speed v 2 . The second and third column illustrate the bias relative to the focal population at the extremes of the interval [0, 1]. Black dots are averages of six replicates, the shaded area represents the 95% confidence interval.

Figure S4 :

 S4 Figure S4: Bias in binary mixes where the speed of the two populations is changed, while maintaining a fixed speed differential. The first column displays the fraction of aggregated particles when each population aggregates separately at the same total packing fraction ρ = 0.07 (other parameters were set as indicated in table S1). Red dots indicate the speed v 1 of the slower, focal population, blue dots the speed v 2 , chosen such that the speed difference ∆v = 4. The second and third column illustrate the bias relative to the focal, slower population at the extremes of the interval [0, 1]. Black dots are averages of six replicates, the shaded area represents the 95% confidence interval.

Figure S5 :

 S5 Figure S5: Fraction of aggregated particles at the end of the simulation as a function of particles adhesive strength and particle speed. The packing fraction was fixed to 0.07 and the remaining parameters of the model were set as indicated in tableS1. The white dotted line represents a fit of the points corresponding to the highest fraction of aggregated particles with Eq. 5.

  signals could be approximated by choosing larger radii, but that would require making additional assumption on signal take-up by neighbouring cells. The intensity of angular noise η and the relaxation time of interaction τ were set such that particles clusters do not display collective motion, similar to what can be observed during early aggregation of Dictyostelium discoideum. Indeed, a low level of angular noise and/or a low relaxation time favor collective motion of particles as interactions contribute to aligning the self-propelled velocity and nematic order emerges. Moreover, we chose a relatively high value of η so that, in the absence of particles interaction, the trajectory of a particle resembles that of a starved D. discoideum cells in the absence of interaction with other cells, i.e. at low cell density (See Fig.1C and Dfor a comparison). Other parameters were set to match the order of magnitude of observations of D. discoideum under standard laboratory conditions: cell size R exp eq = 11µm (De Palo et al., 2017), cell speed v exp = 11µm/min (McCann et al., 2010) and duration of Dictyostelium aggregation t exp f = 12h.

  Variation in particle speed gives rise to different finite-time aggregation patterns in homogeneous populations, leading to the production of non-aggregated cells through two qualitatively distinct routes: delayed aggregation and evaporation. It is thus natural to expect that mixing particles of different speed would reflect into differential rates of aggregation. It is however not obvious how the emergent aggregation patterns, which result from individual interactions, are related to the collective behaviour of the two particle types in isolation. In order to explore this question, we run simulations of binary mixes of particles with two different constant speeds v 1 and v 2 .

Table S1 :

 S1 Model parameters values

	Parameter	value
	N	10000
	t f	100
	dt	0.01
	v 0	8
	µ	1
	η	10
	τ	5
	v 0	8
	F adh	7
	Frep	40
	Req	1.1
	R 0	1.6
	L	364
	ρ = N π(Req/2) 2 /L 2	0.07

Acknowledgments

The authors are grateful to Gabriel Peyré for his help in optimizing the Python program. The authors also thank Thibaut Arnoulx de Pirey Saint Alby for helpful discussions on the link between the adhesion strength and the rate of evaporation and Nirbhay Patil for discussions on active matter physics. This work has received support under the program "Investissements d'Avenir" launched by the French Government and implemented by ANR with the references ANR-10-LABX-54 MEMOLIFE and ANR-10-IDEX-0001-02 PSL* Université Paris, Q-life ANR-17-CONV-6150005, and the project ANR-19-CE45-0002 'ADHeC' PSL research University. L. Brunnet thanks the Brazilian agencies CNPq and Capes for their support.

Supplementary Figures

Supplementary Information

Parameters choice

The fixed parameters of the model were set as indicated in table S1. As in