Sociogenesis in Unbounded Space: Modelling Self-Organised Cohesive Collective Motion

Maintaining cohesion between randomly moving agents in unbounded space is an essential functionality for many real-world applications requiring distributed multi-agent systems. We develop a bio-inspired collective movement model in 1D unbounded space to ensure such functionality. Using an internal agent belief to estimate the mesoscopic state of the system, agent motion is coupled to a dynamically self-generated social ranking variable. This coupling between social information and individual movement is exploited to induce spatial self-sorting and produces an adaptive, group-relative coordinate system that stabilises random motion in unbounded space. We investigate the state-space of the model in terms of its key control parameters and find two separate regimes for the system to attain dynamical cohesive states, including a Partial Sensing regime in which the system self-selects nearest-neighbour distances so as to ensure a near-constant mean number of sensed neighbours. Overall, our approach constitutes a novel theoretical development in models of collective movement, as it considers agents who make decisions based on internal representations of their social environment that explicitly take into account spatial variation in a dynamic internal variable.

Here we discuss details relating to simulations performed to obtain the results presented in the main text.All simulations were performed using a codebase written in the C++ programming language.An asynchronous time updating scheme is used in simulations to avoid parity issues arising from symmetries in agent behaviour.For example, when a synchronous updating scheme is used, two neighbouring agents may exchange positions without ever meeting on the same lattice site.This situation is avoided when using an asynchronous updating scheme.A flow diagram of the timing scheme is shown in Fig. 1.At each discrete time-step, t, the order in which the N agents move is randomised.When an agent, k, is selected to move, they firstly check if there are any other neighbouring agents on the same lattice site, x k .If neighbours are detected, one neighbour, j, is selected at random as an interaction partner.An interaction then takes place between k and j with probability P kj (S k , S j ).In simulation results presented in the main text, agents do not interact if they already had an interaction in the current time-step, or if the currently encountered neighbour was their most recent interaction partner.Note that these conditions mainly influence the system by slowing down the convergence rate of social ranks, and the same stability regimes form without them.The social ranks of the winner and loser are then modulated up and down by δ + and δ − respectively.The initially selected agent, k, then performs a regression measurement with probability r.A high value for r enables agents to be more responsive to socio-spatial dynamics of their neighbours, but incurs a larger computational load for individual agents, and when performing simulations to study the system.On the other hand, when the measurement rate is too low, the rate of information updating is slow, resulting in agents being less responsive to system dynamics.Finally, the agent, k, updates its potential center, x * k , and takes a step according to its individual potential, V k (x; x * k , L).This completes the move of agent k, and the next agent in the randomised list then moves.

Agents on same lattice site as ?
k

Yes
Increase timestep Flow diagram for the timing scheme of simulations.The first time-step begins at the top-left corner of the flow diagram, and proceeds by following the direction of the arrows.The simulation ends when t reaches the pre-specified maximum number of time-step iterations.

A. Polynomial Regression
Agents compute their belief, S k (x), using polynomial regression, a linear regression method that assumes S k (x) to have the form of a polynomial of order u, where ξ is assumed to be a Gaussian white noise with variance σ 2 .Such a belief may also be represented by the coefficient vector β k = (β 0 , . . ., β u−1 ) T .When an agent, k, senses n neighbours, k obtains a data-set consisting of n observations {x j , S j } n j=1 .Linear regression then consists of solving a set of n simultaneous equations, corresponding to Eq. (S1), for the coefficients in β k .Expressing the simultaneous equations in matrix form, this is equivalent to solving where and for simplicity, it is assumed that each of the observations is subject to the same external process noise ξ j = ξ, so that the problem is homoscedastic.The ordinary least squares solution for the coefficient vector, β k , is found by minimising the residual sum of squares (RSS) between observations S j and belief predictions S k (x j ) whose minimum is found through differentiation with respect to the coefficient vector β k , and solving the resulting set of simultaneous equations.The ordinary least squares solutions is then given by [1]

B. Weighted Polynomial Regression
In the model presented in the main text, agents effectively perform local regression [1] to collectively estimate the social field, S(x, t).Weighted polynomial regression is used in order for these estimates to be smooth, by scaling the contribution of each sensed neighbour, j, to the regression by a distance dependent weight factor w j .The weight is defined according to the distance from each neighbour, j, located at x j to the location of the measuring agent, x k , through a smoothing kernel, w j = K h (x k , x j ), defined below.This smoothing kernel makes the agent sensing more local by reducing the contribution of data associated to agents who are further away from the agent performing the measurement.Furthermore, the smoothing kernel, K h , can be considered as a model of agent perception, by defining the extent and strength of agent sensory capabilities.
The smoothing kernel, K h (x k , x j ), has a general functional form [1] where different choices of G define different kernels, and h is the width of the smoothing neighbourhood, which in the present case is the sensory radius of the agent, k.Here, the tri-cube kernel [1] is chosen as it offers compact support, which is desirable for computational convenience.If the support is not compact, agents are assumed to possess an infinite sensory range, which takes into account all other agents in the system for each computation.Weighted regression computes S k (x) by minimising the residual sum of weighted squares between observations, S j , and belief predictions, S k (x j ), given by with normalised weights W jj = w j / k w k .The weighted least squares solution is found by the same method as in the previous section [1].In matrix form, this corresponds to computing the coefficients using the equation where W is a diagonal matrix with elements W jj = w j / k w k , so that j W jj = 1.

C. Model Fidelity
The model fidelity, Φ k , associated to a belief, S k (x), is used by agents as a local measure of fitting accuracy.The quality of the fit is evaluated by considering the amount of information in the sensed data that is explained by the computed belief, by comparing the unweighted residual sum of squares (RSS) to the unweighted total sum of weighted squares (TSS) that would result from a zero-order polynomial regression, S k (x) = β 0 , that is equivalent to an average of S j , where S is the average of all S j in the fitting domain.Supposing the agent, k, located at x k senses n neighbours, the model fidelity is given by which results in a value Φ k ∈ (−∞, 1], with Φ k → 1 indicating a good fit, and Φ k → −∞ indicating a worse fit.Mathematically, the model fidelity is equivalent to computing an unweighted coefficient of determination for a model fitted by weighted regression.Weighted regression model fitting is performed by giving a higher weight to data associated with neighbours that are closer to the sensing agent, and so the estimated model is "constrained" to fit these data points more closely.Hence, by construction, the weighted model may fit the data associated to far-away neighbours more poorly.The model fidelity, Φ k , is evaluated by weighing all data equally, irrespective of distance from the sensing agent.Hence, unlike the coefficient of determination, it is possible that Φ k < 0 when the states of far-away agents diverge greatly from the estimated weighted regression model.Restricting agents from responding to estimates with Φ k < 0 by using a measurement waiting time t c prevents agents from inferring incorrect motion decisions from their weighted regression model when it disagrees greatly with the states of far-away neighbours.

D. Problem Conditioning and Data Quality
For a regression measurement to be performed, the agent must have at least 3 neighbours within sensing distance, h.This is due to the fact that a second-order polynomial regression has 3 coefficients and so requires at least 3 datapoints to fit.Given that a regression is performed, the weighted least squares regression solution in equation (S9) requires a square symmetric moment matrix, X T W X, to be inverted.For this matrix to be invertible, the problem must be well conditioned, as poor conditioning of X T W X may lead to numerical instabilities when it is inverted [2].The reciprocal condition number is a measure of how fluctuations in the input data change the resulting regression solution.This number depends on aspects of the moment matrix such as its rank and scale, with a low reciprocal condition number signifying poor conditioning.In simulations performed, a reciprocal conditioning number threshold, ϵ c = 10 −6 , is defined, below which the moment matrix, X T W X, is considered to be too poorly conditioned to invert.

E. Data Standardisation for Regression
Feature re-scaling [3] improves the regression problem conditioning by mapping the sensed values of x j to a standardised scale.This allows for a single reciprocal conditioning number threshold, ϵ c , to be used, since the scale of the input data has a direct effect on the reciprocal conditioning number.The regression method used assumes sensed data to have a normal distribution [1], and so standardisation is used in order to map the data, x j → x ′ j , so that it has a standard normal distribution, x ′ j ∼ N (0, 1), according to where x and σ x are the mean and standard deviation of the input data, x j .
Due to the linearity of the regression, any re-scaling of the regression coefficients β k → β ′ k that results from rescaling of data x j → x ′ j is reversible after computation of β ′ k .In the case when S k (x) = β 0 +β 1 x+β 2 x 2 , the coefficients for the original regression problem, β k = (β 0 , β 1 , β 2 ), can be retrieved from the re-scaled coefficients, III. NAVIGATION PROCEDURE The potential center locations, x * k (t), are found at each time-step, t, through the solution to the minimisation problem which can be considered graphically as finding the points of intersection between the curve defined by a belief, S k (x), and the horizontal line at , the solutions to Eq. (S16) are given by which yield zero, one or two real solutions.In case two solutions are found, the solution that is closest to the previous potential center location, x * k (t − 1), is chosen, If there are no real solutions, the location, x * k , that minimises Eq. (S16) is the maximum of the belief, S k (x), given by If β 2 = 0, meaning that the agent had performed a first-order polynomial regression, the solution to Eq. (S16) is the location, x * k (t), which solves which is given by Finally, a threshold, ϵ β = 10 −4 , is used, such that the potential center, x * k (t), is only updated when |β 1 | > ϵ β .This is done in order to mitigate large fluctuations in computed potential center locations resulting from near-constant beliefs, S k (x), i.e. when β 2 = 0 and β 1 ≈ 0. This may occur when all agents have similar social ranks, S k , for example in early stages of the formation, or near minima of S(x, t) when agents fit a first-order polynomial regression.In both cases, the low magnitude gradient leads to potential centers, x * k , that can be far away from the current location of the agent, removing them from the group.This parameter requires tuning, since too high a value leads to unresponsive agents, whereas too low a value leads to agents that are too responsive to low magnitude gradients.

IV. INTEGRATED SIGMOID SHAPE FUNCTION
The integrated sigmoid function with respect to x, where K is obtained as a constant of integration.It can be shown that the integrated sigmoid in Eq. (S22) interpolates between an absolute value function (α ≈ 0) and a quadratic function (α → ∞).Using L'Hospital's Rule as α → ∞, When x > 0, the first term tends to zero.When x < 0, the first term tends to 2mx.Hence, When α ≈ 0, F (x; α, m, K) may be expanded around α = 0 using the Taylor series ln(1 + e y ) = ln(2) + where y = −αx.Simplifying the resulting expression, we obtain the quadratic function for α ≈ 0. However, α alone is not enough to classify the different regimes due to the additional dependence of F (x; α, m, K) on m and K, as shown in Fig. 2. Namely, values may be chosen for m and K which counteract the influence that α has in its asymptotic regimes, for example sending m → ∞ or K → ∞ in Eq. (S22).For this reason, the variable J = αK/m − ln(4) is introduced, which captures this multiple asymptotic dependence.This form of J is chosen because it also determines the number of real roots of Eq. (S22), and so determines the region of validity for the fit.
The roots, x * ± , of the integrated sigmoid are found by solving F (x; α, m, K) = 0. Letting z = αx in Eq. (S22), this reduces to Re-arranging for z, exponentiating, and taking a square root both sides results in The inequality in Eq. (S40) is true if the inequality in Eq. (S42) is true.Hence, letting J = Kα/m − ln(4), this defines a lower bound such that the integrated sigmoid has real roots only when J ≥ 0. The value J > 0 corresponds to two real solutions, whereas J = 0 corresponds to one real root where the function maximum touches the x−axis, and the rest of the function is negative.In this latter case, the integrated sigmoid no longer defines a valid fit for the configuration of agents.4) enables us to distinguish the different macroscopic configurations between the different regimes.In both cases, the unstable regime has not been considered, because here agents in the system expand away from each other to such an extent that the integrated sigmoid no longer provides an informative fit.Data corresponds to the results presented in Fig. 3 of the main text.

V. MEASUREMENT WAITING TIME EXPERIMENTS FROM VARYING INITIAL CONDITIONS
To test the efficacy of the two-stage dynamics, we investigate the ability of the system to converge to a concave annular state from a variety of initial conditions when t c = 0 and when t c = 3 × 10 6 .The system is considered to have converged at time t * when the socio-spatial correlation first passes below a pre-defined threshold value, C(t * ) < C * = −0.97.The 'W'-initial condition, shown in the inset of Fig. 3(a), is defined by an initial social rank profile where m = 1, . . ., N is the ranked initial position of each agent from left to right and z 1 = 0, . . ., N/2 interpolates between initial convex (z 1 = 0) and concave (z 1 = N/2) annular configurations.We find that for both conditions of t c , as z 1 is decreased from z 1 = N/2, there exists a threshold value, z * 1 , below which the system no longer converges to a concave annular state.This value of z * 1 is lower for t c = 3 × 10 6 than for t c = 0, indicating that the two-stage dynamics increases the likelihood of reaching the desired goal state.This is illustrated in the insets of Fig. 3(a), which show the system time evolution with and without t c for a selected value of z 1 = 30.When t c = 0, the system remains frozen in the initially defined state and the presence of the radially increasing (convex) edges, facing away from the group center, is maintained.Motion biased by beliefs of agents located at the minima of S(x, t) quickly pushes the system apart so that agents cannot interact with or sense their neighbours.The system becomes locked into the initial configuration, which remains for the duration of simulations.The presence of a measurement waiting time, t c > 0, enables agents to have more interactions in the minima of the configuration.This differentiation in social ranks leads to a destabilisation of the radially increasing edges, which re-configure themselves into radially decreasing, concave edges.
The 'W'-initial condition experiments suggest that the existence of concave, radially decreasing edges at early times is a precursor for the system to reach the desired cohesive pattern at long times.To confirm this hypothesis, we define the concave edge initial condition by initialising the social ranks with S C 0 (z 2 , m) = ∆S max (N/2 − |N/2 − m|) when the ranked position m = 1, . . ., N is at least z 2 agents away from the center (m < z 2 or m > N − z 2 ), and otherwise the agent social rank is randomly sampled from a uniform integer distribution S C 0 (z 2 , k) ∼ U (0, N ∆S max /2), where S k = N ∆S max /2 is the theoretical maximum social rank in the concave annular state.As shown in Fig. 3(b), here a threshold value, z 2 = z * 2 , also exists, below which the system does not converge to an annular state.The value of z * 2 is decreased when t c > 0, illustrating that the two-stage dynamics allows the system to reach the goal state from a larger perturbation.Insets correspond to z2 = 3.In this case, the configurations obtained for tc = 0 are concave but not radially symmetric.In both cases, the inclusion of tc > 0 increases the range of initial condition states that converge to a concave annular state, defined by a socio-spatial correlation C(t) < −0.97.Simulation parameters are the same as those in Fig. 2 of the main text, with k h = 0.5 and kL = 1, and results are averaged over 50 Monte Carlo realisations.

A. Social Convergence
The systems exhibits convergence to a concave annular state in the stable Full Sensing and Partial Sensing regimes, as shown in Fig. 4(a).It can be seen that parameter values with lower magnitude socio-spatial correlation are found in the Unstable regime, where the sensed proportion of agents is below n s (t) < 0.06, which in simulations corresponds to approximately 6 sensed agents.Below this value, the collective agent motion resulting from computed beliefs, S k (x), becomes unstable.If such instabilities lead small groups of agents to move away from the group, the group center of mass is shifted, thereby reducing the magnitude of the socio-spatial correlation.

FIG. 2 .
FIG. 2. (a)The shape regimes of the integrated sigmoid in Eq. (S22) cannot be classified using α alone.(b) Using J = αK/m − ln(4) enables us to distinguish the different macroscopic configurations between the different regimes.In both cases, the unstable regime has not been considered, because here agents in the system expand away from each other to such an extent that the integrated sigmoid no longer provides an informative fit.Data corresponds to the results presented in Fig.3of the main text.

FIG. 3 .
FIG. 3. (Colour online) Convergence times, t * , from different initial conditions with and without a measurement waiting time, tc.(a) 'W'-initial condition, S W 0 (z1, m).Insets correspond to z1 = 30.(b) Concave edge initial condition, S C 0 (z2, m).Insets correspond to z2 = 3.In this case, the configurations obtained for tc = 0 are concave but not radially symmetric.In both cases, the inclusion of tc > 0 increases the range of initial condition states that converge to a concave annular state, defined by a socio-spatial correlation C(t) < −0.97.Simulation parameters are the same as those in Fig.2of the main text, with k h = 0.5 and kL = 1, and results are averaged over 50 Monte Carlo realisations.

4 FIG. 4 .
FIG. 4. (a)Socio-spatial correlation, C(t), in the (kL, k h )-state-space plotted after t = 10 9 steps.Lower magnitude values of C(t) are obtained in the Unstable regime, where unstable expansion leads some agents to move away from the group.This shifts the group center of thereby the socio-spatial (b) Scatterplot of SSC plotted against the average proportion of sensed neighbours ns(t).Points are coloured by the corresponding value of J and their shape is determined by the mean nearest-neighbour distance ∆(t).Colourless points are those with ns(t) < 0.06 and correspond to the Unstable regime.Data corresponds to the results presented in Fig.3of the main text.

FIG. 5 .
FIG. 5. (a)Fitted power-law coefficient, |γν |, for the mean-squared-displacement, ν(t), at time t = 10 9 .The absolute value is taken so that data can be plotted on a logarithmic scale when γν is close to zero and negative.Higher power-law coefficient values are found on the cusp between the Full and Partial Sensing regimes, where the system exhibits two possible expansion states and has not yet reached equilibrium.The power-law coefficient was fitted between t = 6.45 × 10 8 and t = 9.95 × 10 8 .(b) Corresponding MSD curves plotted on a log-log scale for each parameter-pair kL and k h .On each vertical axis, the MSD, normalised by the initial value is plotted, while time, t, is plotted on the horizontal axes.Data corresponds to the results presented in Fig.3of the main text.

FIG. 6 .
FIG. 6. (a)Fitted power-law coefficient, |γ∆|, of the mean nearest-neighbour distance, ∆(t), at time t = 10 9 .The absolute value is taken so that data can be plotted on a logarithmic scale.As in Fig.5, higher magnitude values of |γ∆| are found on the cusp between the Full and Partial Sensing regimes.Power-law coefficient fitting times are the same as in Fig.5.The colourbar has a lower limit of |γ∆| > 10 −3 to aid visualisation, as values in the Unstable regime reach |γ∆| = 10 −15 .(b) Normalised mean nearest-neighbour distance for each parameter-pair, kL and k h .Vertical axes show ∆(t) normalised by the initial value, while time, t, is plotted on the horizontal axes (log-log scale).Data corresponds to the results presented in Fig.3of the main text.