Emergent order in a continuous state adaptive network model of living systems

Swarming is a collective phenomenon which is encountered in many biological systems, such as schools of fish, swarms of locusts and flocks of birds. Some swarming phenomena in biological systems can be guided by external cues, or a leading individual, but a large variety of biological systems display spontaneously emerging, self-organized group behavior [1]. Even unicellular organisms, such as the bacterium E. coli and the amoeba D. discoideum, have been observed to display such collective behavior [2].

Models of swarm formation can be broadly put into two categories. The first group we call mechanistic models. Their microscopic interactions are solely based on first principles, and they do not encode the tendency to order explicitly. This works well if all microscopic forces that determine the dynamics are known. Such models have been successfully used to describe the dynamics of self propelled particles of different shapes that exhibit exclusion interactions [3, 4, 5]. The second group, that we call heuristic models, are well exemplified by opinion formation in social groups. We know that people can convince others to join their cause. However, we do not yet understand the details of how one individual can cause another to form an opinion, thus we have a need for a heuristic rule that postulates how two individuals can form a consensus. Such heuristic treatment is in general needed when it is clear that there is a microscopic tendency to align, but we do not know the details of interactions; in context of opinion dynamics, such a model was first put forward by Vicsek et al. [6].

Giving up the connection to the mechanistic first principles opens a door to an interesting opportunity. One can recognize that the heuristic of be more like your neighbor is now the central aspect of the evolution of the system, instead of spatial dynamics of the particles. Thus we are shifting our attention from spatial dynamics to the topology of interactions, a problem for which the language of complex networks is uniquely well suited. In this framework the system is represented by a dynamic network. The agents, corresponding to the nodes of the network, can occupy different states, representing the direction of movement, and the changing topology of the network encodes the temporal dynamics of the interactions. An example of such a system can be seen in fig. 1. One such dynamic network model, for a system with a discrete state space, has been proposed by Chen, Huepe and Gross [7]. This model predicts a first or second order phase transition from a disordered to an ordered state, depending on the number of possible internal states. The first order phase transition has previously been observed in agent based models that follow the spatial dynamics of individuals, and has been confirmed both in experiments with active matter and observations of living systems [8, 9, 10]. The character of the transition can also depend on finite system size effects [11, 12] and on subtle changes in the way the noise is incorporated [13].

Earlier network models [14, 15, 7] have focused on a discrete state space. All the aforementioned biological systems, however, have a continuous state space. Our intuition from equilibrium statistical mechanics tells us that the nature of the state space can matter a lot [16, 17, 18]. It has been demonstrated that in nonequilibrium systems, details affect the type of phase transition. Therefore, in this paper we investigate the swarming transition in a continuous state space.

In this letter we present the numerical solutions of our model accompanied with analytic expressions for the phase space boundary, that we found through stability analysis of the disordered phase. Our results clearly show that continuous state adaptive network models are able to describe spontaneous emergence of order in active systems. We find that the balance of time scales in the model determines the type of phase transition. Our model has two main time scales, the time scale of information propagation through the network, and the time scale of network reshaping. When the latter is significantly faster than the former, we get a mean field like situation, where the network is updated so fast that every node has contact with a random selection of nodes, replacing the link structure by an average connectivity number. The phase transition in this case becomes second order, while in the case that the time scales are of comparable magnitude, we find a first order transition. A meta phase transition between the first and second order regimes, similar to this one has been reported before in the context of active adaptive systems [12].

Model – We model a system of self-propelled particles with a constant speed and changing direction in two dimensions. Each particle corresponds to a node in a network, with an internal state that represents the direction of movement. We represent this internal state by an angle $\theta\in(-\pi,\pi]$. Nodes may be connected by links, indicating mutual awareness. We will refer to individuals connected by a link as neighbors.

In our system not only the states of individual agents can evolve dynamically, but also the relationship of the agents with their surroundings. Analogously to [7], we distinguish four types of dynamics, given by four rules:
Rule 1 Individuals spontaneously change their heading direction to another uniformly chosen direction with rate $w_{0}$.
Rule 2 Individuals adopt to the average direction of two neighbors with rate $w_{2}$.
Rule 3 Arbitrarily chosen not-neighboring individuals become neighbors with a coupling rate $c$.
Rule 4 Arbitrarily chosen neighbors loose mutual awareness with a decoupling rate $d$.

An illustration of the model is given in fig. 2. The first two rules are comparable to the rules in the Vicsek model [6]; rule 1 is similar to the noise that is added to each particle’s updated heading direction, whereas rule 2 corresponds to the tendency of individuals to align with their neighbors within a certain radius. The radius is modeled in the network with the links, since we do not keep track of the physical position of individuals in space. The interactions described by rule 3 and rule 4 are needed since non-neighboring individuals, moving in different directions, might become aware of each other and conversely individuals which are neighboring, but head in different directions, may loose mutual awareness; therefore, their link must be created or removed respectively. Since we want to look at the behavior of groups comprised of agents of the same species, we assume that the rates are global.

We define the state distribution function $f(\theta;t)$ which represents the fraction of the agents in the state $\theta\in(-\pi,\pi]$ at time $t$. Thus the integral of $f(\theta,t)$ over the whole domain is normalized at any time. The cartoon in fig. 1 is an illustration of such a state density that starts out as a random, noisy distribution over the states and evolves into a state corresponding to a collective motion in one predominant direction. Moreover, we define the link distribution function $l_{2}(\theta,\theta^{\prime};t)$ as the density of neighboring individuals, one of which is heading in direction $\theta$ while the other one has direction $\theta^{\prime}$, at time $t$. Note that in contrast to $f(\theta;t)$, the link density $l_{2}(\theta,\theta^{\prime};t)$ is not normalized since the overall number of links can grow or shrink, making the normalization of $l_{2}$ time dependent. We also define higher order interaction terms like the path-like three node and star-like four node subgraph densities $l_{3}(\theta;\theta^{\prime},\theta^{\prime\prime};t)$ and $l_{*}(\theta;\theta^{\prime},\theta^{\prime\prime},\theta^{\prime\prime\prime};t)$ analogously. A visual representation of these subgraphs can be seen in fig. 2. We note that the link density functions obey certain symmetries. In particular, we have $l_{2}(\theta,\theta^{\prime};t)=l_{2}(\theta^{\prime},\theta;t)$ and for higher order terms all non central nodes can be exchanged in any order.

One can derive the master equations for the continuous state model using the four rules, which govern the dynamics of the state and link density distributions. These equations can be interpreted as the extension of the master equations in [7], into the continuous state set.

The first term of eq. 1a scales with $w_{0}$ and models the spontaneous direction changes of nodes. $F^{\mathrm{int}}$ in the second term is a functional that describes three-body interactions that cause a change of the state of an agent. The eq. 1b contains four contributions. $L^{\mathrm{noise}}$ accounts for changes in link density, due to the random changes of states of already linked agents. $L^{\mathrm{int}}$ represents the changes of the link density due to the changes of states of already linked agents, but in this case the changes were induced by interactions between neighbors. The last two terms represent the change in link density due to link creation or deletion and scale with the rates $c$ and $d$ respectively. Both $L^{\mathrm{noise}}$ and $L^{\mathrm{int}}$ integrate to zero over the whole domain [19] i.e. $\int_{-\pi}^{\pi}\mathrm{d}\theta\;\int_{-\pi}^{\pi}\mathrm{d}\theta^{\prime}\ L=0$, thus we can show that the overall link density, that we define as

always evolves to a steady state. Integrating eq. 1b over $\theta$ and $\theta^{\prime}$ yields a simple differential equation $\partial_{t}k(t)=c-dk(t)$ which has the solution

Thus we find for the steady state link density $k_{s}=c/d$. During our analysis we found that the system was not qualitatively affected by changes in the value of $k_{s}$ [19], so from now on, unless explicitly stated otherwise, we will assume that $k_{s}=1$, i.e. $c=d$. We will refer to both rates as the rate of link dynamics. The parameter space can be reduced further by rescaling time with the rate $w_{2}$. Going forward we will set $w_{2}=1$, which implies that time is measured in units of $1/w_{2}$ and all other rates are measured in multiples of $w_{2}$.

In order to solve eq. 1 we need to relate higher order subgraph density functions $l_{3}$ and $l_{*}$ to the link density function $l_{2}$

We call the resulting model the moment closure approximation (MCA) model. The MCA model can be solved directly, and we can obtain a phase space, with ordered and disordered phases. In the case of a high rate of link dynamics, i.e. $d\to\infty$, the time evolution of the link density in eq. 1b leads to the steady state of

Equation 6 shows that for large sampling rates, the probability of two states being linked is proportional to the relative abundances of those states. In this case eq. 1a becomes independent of the link density and reduces to a single differential equation for $f(\theta;t)$. We recognize this as the mean field approximation of the full model. We call the steady state solution of the mean field model $f_{s}(\theta)$, which is disordered if all states occur with same frequency i.e. $f_{s}(\theta)=\frac{1}{2\pi}$. We call the solution ordered if one state is more abundant than all others, see fig. 1. To properly quantify this notion of order we introduce an order parameter $a_{1}$, which is the amplitude of the first symmetric Fourier mode in the series expansion of $f_{s}(\theta)$, given by

In eq. 7 the $a_{n}$ are the mode amplitudes and $\theta_{s}$ is the direction preferred in the steady state. Since $\cos(\theta-\theta_{s})$ is a function with a single peak around $\theta_{s}$, the contribution of the first mode to the series expansion of $f_{s}(\theta)$ gives us a quantitative understanding of how ordered the system is. The atypical normalization of the Fourier series is chosen to have $0\leq a_{1}\leq 1$. Thus $a_{1}=0$ is associated with the disordered and $a_{1}=1$ with the ordered state.

Results – We first solved the model in the mean field approximation, since in this case we only have an equation for $f(\theta;t)$ to solve, with $w_{0}$ as its only control parameter. Using the Fourier series expansion in eq. 7, we can reduce the differential equation for the mean field system to a set of algebraic equations for the Fourier modes in steady state

where $\gamma_{npq}=$
$\int_{-\pi}^{\pi}\mathrm{d}\theta\int_{-\pi}^{\pi}\mathrm{d}\phi\cos\left[p\left(\theta+\frac{\phi}{2}\right)\right]\cos\left[q\left(\theta-\frac{\phi}{2}\right)\right]\cos\left(n\theta\right).$ From eq. 8 we can extract the exact value of the critical noise

and the scaling behavior of the order parameter close to the critical point $a_{1}\propto(w_{0}-w^{\mathrm{cr}}_{0})^{1/2}$. Our analytic approximation of $a_{1}$ is obtained by closing eq. 8. This procedure provides a generic $n$-th order polynomial for $a_{1}$ if closed at that order. Thus the best possible analytic solution is a root of a fourth order polynomial, which is plotted in fig. 3 as a solid black line. This solution gets progressively worse away from the critical point, but close to it our numerical and analytic results are in great agreement. Furthermore, we see that this phase transition is unambiguously second order and has a critical exponent of $1/2$, which was obtained analytically and is exact [19].

For the MCA model we have an additional control parameter $d$, which sets the time scale of the link dynamics. Therefore, we get a critical line

which can be seen in fig. 4, as the red dash-dotted curve. The vertical red dashed line is drawn at the value $w_{0}^{\mathrm{cr}}$, this value is reached when $d$ goes to infinity. This critical line is obtained by linear stability analysis of the disordered state and is exact, just like the mean field model. The disordered state becomes unstable left of the critical line. Our numerical solutions were obtained by solving the system from an ordered initial condition. Thus the numerical solutions approach the phase transition from the ordered side. The existence of numerical results with nonzero order parameter, on the right of the line indicates hysteresis. Unfortunately, the methods that allowed us to obtain analytic values of the order parameter near criticality in the mean field model do not work here. However the numerical results in fig. 4, clearly show a discontinuous change in the value of the order parameter at the critical line. Together with the bistable region, this is a clear sign of a first order phase transition. To exclude artifacts related to slow convergence to steady state, we repeated the runs for $3\times 10^{5}$ and $9\times 10^{5}$ time steps measured in units of $1/w_{2}$. The results remained unchanged [19].

Discussion – The continuous state adaptive network model can describe the emergence of swarming. The key parameter in this transition is the relative strength of the noise compared to the agent interactions. The critical noise value, at which the swarming transition occurs, can be changed by tuning the timescale of the network dynamics. However, if this timescale, given by $1/d$, becomes much smaller than the timescale of the node dynamics, set by $1/w_{2}$, two interesting things happen. First, the region of stability shrinks and the swarming phase transition changes from a first, to a second order transition. The shrinking of the stable region, for fast link dynamics, is surprising because, the limit of this regime is the mean field approximation. Using the intuition from equilibrium statistical mechanics, one would expect the critical noise to be higher in the case of the mean field model, since we expect mean field models to overestimate the tendency to order. However, we get the opposite result. The critical noise is maximal for $d=0$; most likely due to the decoupling of the link dynamics from the state dynamics. This decoupling can slow down the propagation of orientational information through the system. The second striking effect seen in our system is the meta phase transition. While changing the rate of the link dynamics, the swarming transition becomes second order. Further investigation is necessary, but it seems that this effect can be explained by the observation that the type of noise affects the type of phase transition. According to Pimentel et al. [20] intrinsic noise in agent states leads to a second order transition, while extrinsic noise, caused by imperfect sampling of neighbors, leads to a first order phase transition. This fits well with our findings, since the mean field system only has intrinsic noise caused by spontaneous direction change, while the MCA system introduces an additional noise source, through changing topology. We do not know if this meta transition happens at a finite value of $d$ or strictly at $d\rightarrow\infty$, but regardless of the point where the meta transition happens, it constitutes a significant change in the behavior of the system. Biological systems that are described by this model, would therefore not always need nucleation points, such as small co-moving groups, for global swarms to emerge. If such systems could facilitate a quickly changing neighborhood topology, they would be able to smoothly transition into a swarm state.
The new model that we have put forward in this letter demonstrates that the timescale separation between the topology- and state-dynamics of the network has a strong impact on the nature of the swarming transition.