## Abstract

Clustering is a phenomenon that may emerge in multi-agent systems through self-organization: groups arise consisting of agents with similar dynamic behaviour. It is observed in fields ranging from the exact sciences to social and life sciences; consider, for example, swarm behaviour of animals or social insects, the dynamics of opinion formation or the synchronization (which corresponds to cluster formation in the phase space) of coupled oscillators modelling brain or heart cells. We consider a clustering model with a general network structure and saturating interaction functions. We derive both necessary and sufficient conditions for clustering behaviour of the model and we investigate the cluster structure for varying coupling strength. Generically, each cluster asymptotically reaches a (relative) equilibrium state. We discuss the relationship of the model to swarming, and we explain how the model equations naturally arise in a system of interconnected water basins. We also indicate how the model applies to opinion formation dynamics.

## 1. Introduction

The basic principles explaining the emergence of one (or several) clusters in multi-agent systems, as documented in the literature, invariably adopt a simple model for the behaviour of the agents.

Swarming models mostly focus on the collective behaviour and the cohesion of a single cluster maintained by the attraction between the animals (or the alignment of their velocities), in counterbalance with the repulsive interactions and/or the drift induced by random walk behaviour (Vicsek *et al*. 1999; Gazi & Passino 2003; Mogilner *et al*. 2003).

Models for opinion formation often focus on consensus reaching (Hegselmann & Krause 2002), i.e. the process of forming a single cluster of people agreeing unanimously on a particular issue, or the coexistence of only two opposite opinions, usually involving nearest neighbour interactions (Sznajd-Weron & Sznajd 2000; Bordogna & Albano 2007). However, other models have been proposed that produce multiple clusters with different opinions as a general outcome (Deffuant *et al*. 2002).

A particular type of clustering is observed in systems of coupled oscillators, such as aggregations of flashing fireflies and coupled Josephson junction arrays (Strogatz 2003). (The term clustering is then often replaced by *synchronization* or *entrainment*.) One distinguishes between *phase* clustering and *frequency* clustering. The first form (e.g. Golomb *et al*. 1992; Okuda 1993; Zanette & Mikhailov 2004) is associated with networks of oscillators with identical natural frequencies. (The natural frequency of an oscillator characterizes its behaviour in isolation, i.e. without interaction with other oscillators.) Each cluster consists of oscillators with (asymptotically) equal phases, with the number of different clusters depending on the interaction. For non-identical natural frequencies for which the differences are sufficiently small, it may still be possible to distinguish different phase clusters (Tass 1997). Larger differences between the natural frequencies may induce oscillators having different long-term average frequencies, resulting in frequency clustering (Morelli *et al*. 2005; De Smet & Aeyels 2007); each cluster is characterized by the long-term average frequency of its members. For more details on both phenomena and for examples of clustering in chaotic systems, we refer to Manrubia *et al*. (2004).

Due to the complexity and richness of the dynamics of some of these models, analytical results are often hard to come by and are restricted to the existence and local stability properties of some of their solutions (Golomb *et al*. 1992; Tass 1997; Aeyels & Rogge 2004), and exploration of the parameter space is usually carried out by simulations (Grégoire *et al*. 2003; Maistrenko *et al*. 2004).

In previous papers (Aeyels & De Smet 2007, 2008), we introduced a simple model for clustering, corresponding to a system of non-identical attracting agents. It may be considered as a simplification of the Kuramoto (1984) model of coupled oscillators that retains its (frequency) clustering behaviour, while exhibiting an increased potential for developing analytical results. Each agent has an autonomous component—its natural velocity—and attracts other agents by a saturating interaction function. The balance between autonomous behaviour and attraction depends on the value of the coupling strength. For high values of the coupling strength, distances between agents remain bounded, resulting in a single cluster. For lower values, several clusters arise; each cluster is characterized by the common asymptotic velocity of its agents. When the coupling strength equals zero, all agents move at their natural velocities. In contrast to many of the aforementioned models, a thorough analysis is possible. We were able to analytically characterize the cluster structure by a set of *necessary and sufficient* inequalities in the model parameters and we showed that there exists a *unique* cluster structure satisfying these inequalities for a given set of parameters and that the distances between agents from the same cluster approach constant values.

In the last decade, researchers have introduced various network models representing links between agents in an attempt to reproduce some of the characteristics of real-life networks. Based on a few simple rules, one is able to generate networks representing collaborations between film actors or scientists, links in the World Wide Web or social networks, which are called, respectively, small-world networks (Watts & Strogatz 1998), scale-free networks (Barabási & Albert 1999) and community networks (Jin *et al*. 2001). For an overview, we refer to Albert & Barabási (2002) and Newman (2003). An important characteristic of small-world networks and community networks also involves the notion of ‘clustering’, of course with a different interpretation since no dynamics are associated to the nodes.

A question that naturally arises is how the behaviour of dynamical systems—and in particular any emerging clustering behaviour—is affected by the network structure that defines the interactions. Research in this direction has mostly focused on synchronization in networks of coupled dynamical systems (in particular, in systems of coupled oscillators; Hong *et al*. 2002; Li & Chen 2003; Jadbabaie *et al*. 2004; Moreno & Pacheco 2004). The present paper also fits into this research direction.

In this paper, we show that the aforementioned results of Aeyels & De Smet (2007, 2008) are not restricted to the all-to-all coupled case, but can be extended to more general interactions with respect to both the interaction function and, perhaps more importantly, with respect to the underlying network structure. To illustrate the potential of the extended model, we indicate how the model equations naturally arise in a system of interconnected water basins (generalizing an example discussed in Aeyels & De Smet 2008) and we illustrate its relevance to swarming and opinion formation dynamics.

In the following sections, we explain the model, define the concept of clustering behaviour and calculate the asymptotic average cluster velocity. We then present and discuss the results of this paper and provide an outline of the proofs. In §§7–9, we show how the model applies to swarming, compartmental systems and opinion formation, respectively. For the detailed proofs of the analytical results, we refer to the electronic supplementary material.

## 2. The model

We consider the system described by the following differential equations:(2.1)with *A*_{i}, *γ*_{j}>0, *K*≥0 and *N*>1. For all *i*, *j*∈*I*_{N}, the functions *f*_{ij} are non-decreasing, Lipschitz continuous and satisfyfor some *F*_{ij}∈. (It follows that the functions *f*_{ii} (*i*∈*I*_{N}) are odd, and since they are only evaluated in zero, they play no role, but are retained for notational convenience.) The interval [−*F*_{ji},*F*_{ij}] covers the range of the interaction of agent *j* with agent *i*. If both *F*_{ij} and *F*_{ji} are positive, then the interaction between the agents *i* and *j* will be attractive, except when the distance *x*_{j}(*t*)−*x*_{i}(*t*) between the agents *i* and *j* belongs to the interval between zero and . The extent to which each individual agent *j* tends to interact with other agents is represented by the weight *γ*_{j}. The parameter *A*_{i} reflects the sensitivity of the agent *i* to interactions with other agents and *K* is the global coupling strength.

The conditions that *K*≥0 and that the functions *f*_{ij} are non-decreasing are important to be able to show that the clustering behaviour (defined in §3) is independent of the initial condition. If these conditions are not satisfied, some of the results presented in this paper may still partially hold (e.g. we suspect that the system will still exhibit clustering behaviour), but independence of initial conditions and uniqueness of the cluster structure for a given set of parameters cannot be maintained.

Since we can rewrite the system (2.1) in the form(2.2)by setting ( retains the properties of *f*_{ij}), (*i*, *j*∈*I*_{N}) and then omitting the accents, we will consider only the system (2.2) for the mathematical analysis, corresponding to putting *A*_{i}=1, for all *i* in *I*_{N}.

The model in Aeyels & De Smet (2007, 2008) corresponds to an all-to-all interaction structure with and *A*_{i}=1, for all *i* in *I*_{N}, and where the interaction functions *f*_{ij} are all equal to the same function *f*, which is assumed to be odd, and reaches its saturation valuefor some *F*, *d*>0. As shown in the present paper, the results from Aeyels & De Smet (2007, 2008) can be generalized to the system (2.2). The generalization requires a slightly amended definition of the notion of cluster structure, in comparison with Aeyels & De Smet (2007, 2008). The formulation of the conditions determining a cluster structure also needs some adaptation. The proofs of the results tend to be quite technical at places, mainly due to the general nature of the interactions. In turn, however, new applications may now be envisaged, as will be illustrated in §§8 and 9.

## 3. Clustering behaviour

We consider the following definition of *clustering behaviour* of a solution *x* to (2.2) with respect to a cluster structure *G*=(*G*_{1}, …, *G*_{M}), representing an ordered set partition of *I*_{N}.

The distances between the agents in the same cluster remain bounded (i.e. |

*x*_{i}(*t*)−*x*_{j}(*t*)| is bounded for all*i*,*j*∈*G*_{k}, for any*k*∈*I*_{M}, for*t*≥0).For any

*D*>0, there exists a time after which the distances between agents in different clusters are at least*D*.The agents are ordered by their membership to a cluster:

*k*<*l*⇒*x*_{i}(*t*)<*x*_{j}(*t*), ∀*i*∈*G*_{k}, ∀*j*∈*G*_{l}, ∀*t*≥*T*, for some*T*>0.

This definition differs slightly from the one in Aeyels & De Smet (2007, 2008), but the results from Aeyels & De Smet (2007, 2008) remain valid with this modified definition.

In Aeyels & De Smet (2008), we show that there exist inequalities in the parameters of the model, which constitute a *necessary and sufficient* set of conditions for clustering behaviour of *all* solutions of the all-to-all coupled model with equal weights with respect to a given cluster structure (*G*_{1}, …, *G*_{M}). We also show that, for every given set of parameters *b*∈, *K*∈ and , there exists a unique ordered set partition (*G*_{1}, …, *G*_{M}) of clusters satisfying these conditions, and that the internal behaviour of the clusters asymptotically reaches an equilibrium situation.

We will generalize these results to the model (2.2) (although, for convenience, some non-generic cases are not included in the generalization).

## 4. Asymptotic average cluster velocity

Set(Note that lim_{D→+∞}*φ*(*D*)=0.)

Consider a non-empty set *G*_{0}⊂*I*_{N}. For a vector *w*∈, we define as the weighted average of *w*_{i} over *G*_{0}, with weighting factors *γ*_{i},Owing to the properties of the interaction functions *f*_{ij}, all internal interactions (i.e. interactions between agents in *G*_{0}) will be cancelled in the expression for . Pick *D*>0 and assume that, at some time instance *t*_{0}, the agents in *G*_{0} are separated by at least a distance *D* from all other agents (i.e. |*x*_{i}(*t*_{0})−*x*_{j}(*t*_{0})|≥*D* whenever *i*∈*G*_{0} and *j*∉*G*_{0}). Then the interaction of agent *i* in *G*_{0} with an agent *j* not in *G*_{0} will deviate at most max(*F*_{ij}−*f*_{ij}(*D*), *F*_{ji}−*f*_{ji}(*D*)) from its saturation value. Assume, furthermore, that at *t*_{0} each agent not belonging to *G*_{0} has an *x*_{i}(*t*_{0}) value either smaller or larger than *all x*_{i}(*t*_{0}) values of agents in *G*_{0}. Denote by (or ) the set of agents with *x*_{i}(*t*_{0}) values smaller (or larger) than the *x*_{i}(*t*_{0}) values of the agents in *G*_{0}. Then for any *i*∈*G*_{0}, it follows thatand thus(4.1)where we define the function asfor all , *G*_{0}, ⊂*I*_{N} with *G*_{0} non-empty.

For an ordered set partition *G*=(*G*_{1}, …, *G*_{M}) of *I*_{N}, let denote and let .

If the solution *x* exhibits clustering behaviour with respect to *G*=(*G*_{1}, …, *G*_{M}), then the application of (4.1) with *G*_{0}=*G*_{k}, and (for some *k*∈*I*_{M}), and considering the limit *D*→+∞ (implying *t*_{0}→+∞ owing to the clustering behaviour), results inIn theorem 5.5, it is shown (for the generic case) that this is then also the asymptotic velocity for each member of *G*_{k}.

## 5. Results

We present our main results. The proof is outlined in §6; detailed proofs are given in the electronic supplementary material.

### (a) Necessary and sufficient conditions for clustering behaviour

For any set *S*, let () denote the set of all ordered partitions of *S* in two subsets, i.e.

*Let G*=(*G*_{1}, …, *G*_{M}) *be an ordered set partition of I*_{N}. *The following set of inequalities is* sufficient *for clustering behaviour with respect to* *G of* all *solutions of the system* (*2.2*):(5.1a)*and*(5.1b)*The following set of inequalities is* necessary *for clustering behaviour with respect to G of* all *solutions of the system* (*2.2*):(5.2a)*and*(5.2b)

If the saturation values *F*_{ij} are attained, i.e. if *f*_{ij}(*ξ*)=*F*_{ij} for *ξ* sufficiently large, for all *i*, *j*∈*I*_{N}, then the conditions (5.1*a*) and (5.2*b*) are *necessary and sufficient* for clustering behaviour with respect to *G* of all solutions of the system (2.2). This is easily shown by adaptation of the proof of theorem 5.1, given in the electronic supplementary material; it encompasses a similar result in Aeyels & De Smet (2007, 2008) for the case of all-to-all interactions. The conditions (5.1*a*) and (5.2*b*) are not necessary and not sufficient in case saturation is not attained, as follows from the example below.

An example, showing that, in general, (5.1*a*) and (5.1*b*) are not necessary and (5.2*a*) and (5.2*b*) are not sufficient, is constructed as follows. Consider the system (2.2) with two agents, with *K*=1, *b*_{1}=−1 and *b*_{2}=1, *γ*_{1}=*γ*_{2}=1 and *F*_{12}=*F*_{21}=1. Then (5.1*a*) and (5.1*b*) are not satisfied for any cluster structure, while (5.2*a*) and (5.2*b*) are satisfied for both the cluster structures ({1, 2}) and ({1}, {2}), sinceBoth ({1, 2}) and ({1}, {2}) are possible cluster structures, but which cluster structure characterizes the behaviour of (solutions of) the dynamical system cannot directly be derived from theorem 5.1. However, if the interaction function *f*_{12} reaches its saturation value *F*_{12}, then the above remark guarantees the emergence of the cluster structure ({1,2}), as it satisfies the conditions (5.1*a*) and (5.2*b*). Indeed, subtracting the system equations for agents 1 and 2 results inshowing that the system will exhibit clustering behaviour, with agents 1 and 2 belonging to the same cluster if and only if *f*_{12} attains its saturation value *F*_{12}=1 for finite values of its argument. The other cluster structure ({1},{2}) emerges when saturation is not reached.

### (b) Clustering behaviour with varying coupling strength

*For every b*∈ *and every matrix F*∈()^{N×N} *with F*_{ij}+*F*_{ji}≥0 *for all i*, *j*∈*I*_{N}, *there exists a partition of* *in a finite number of intervals*, *such that each interval corresponds to a unique ordered set partition G of I*_{N}, *for which* (*5.1a*) *and* (*5.1b*) *hold for all K in the interior of this interval*.

When (5.1*a*) and (5.1*b*) are satisfied for values of *K* in the interior of some interval, it follows that (5.2*a*) and (5.2*b*) are satisfied for *K* in the closure of this interval. At the endpoints, more than one cluster structure will satisfy (5.2*a*) and (5.2*b*), and no cluster structures will satisfy (5.1*a*) and (5.1*b*), as is illustrated by example 5.3.

Considering the relevance of theorem 5.4 to system (2.2), the conditions *F*_{ij}+*F*_{ji}≥0 are required for the existence of non-decreasing interaction functions *f*_{ij} with saturation values −*F*_{ji} and *F*_{ij}. In the proof of theorem 5.4, the uniqueness of the clustering behaviour with respect to the initial conditions of solutions to (2.2) guarantees the uniqueness of the cluster structure *G* for a given value of the coupling strength. Without these conditions, the result remains valid (as follows from the proof), except for the uniqueness of the cluster structure.

If all entries of *F* are non-negative and *F* is irreducible (i.e. for any partition (*G*_{1}, *G*_{2})∈(*I*_{N}) there exist *i*∈*G*_{1} and *j*∈*G*_{2} with *F*_{ij}≠0), then for large *K* values one can check that the cluster structure will be equal to a single set containing the entire population, i.e. *G*=(*I*_{N}).

As an illustration of the varying cluster structure for different ranges of *K*, figure 1 shows the long-term average velocities (which are of course equal for agents from the same cluster, and thus the vector *v* reflects the emerging cluster structure) of the agents for varying coupling strength *K*, with , for all *i* in *I*_{6},From this figure, the intervals for *K* and the corresponding cluster structures can be inferred immediately. Note that the clusters may split up with increasing coupling strength; a similar phenomenon can be observed in the all-to-all coupled Kuramoto model (De Smet & Aeyels 2007).

### (c) Internal behaviour of the clusters

To investigate the internal behaviour of a cluster, we will exclude interaction functions that are constant in an infinite number of disjunct intervals of non-zero length, i.e. we only consider functions *f*_{ij} for which the set(with |*S*| denoting the number of elements of a set *S*) is finite. Furthermore, not all, but *almost all* values of *b*∈ are allowed; for convenience, the formulation of the exact conditions under which theorem 5.5 is true is postponed to the proof.

*Consider the system* (*2.2*) *with interaction functions f*_{ij} *having a finite number of function values* *for which* *is an interval of non-zero length*. *For almost all b*∈, *the following is true*.

*If x is a solution with cluster structure G*=(*G*_{1}, …, *G*_{M}), *then for each k*∈*I*_{M}, *if i*,*j*∈*G*_{k},

The conditions on *b* are related to the conditions on the functions *f*_{ij}, e.g. if all functions *f*_{ij} are *increasing* on , then the result of theorem 5.5 holds for *all b*∈, as will be clear from the proof.

The second result of theorem 5.5 (stating that the velocities of the agents approach a constant value) is valid for *all b*∈, and for *all* interaction functions satisfying the conditions introduced in §2. (For a proof, we refer to the future work of De Smet & Aeyels in preparation *a*.)

To show that the first result of the theorem 5.5 does not always hold, we present an example for which the distance between two agents from the same cluster does not converge.

We will construct an example with *N*=5 and *K*=1, where *x*_{4}(*t*)−*x*_{3}(*t*) and *x*_{3}(*t*)−*x*_{2}(*t*) do not converge, although agents 2, 3 and 4 belong to the same cluster. All weighting factors *γ*_{i} are set equal to 1. The vector *b* and the matrices *F* and *d* (where ) are chosen asThe resulting cluster structure is *G*=({1}, {2, 3, 4}, {5}), with asymptotic velocities, respectively, −1, 0 and 1, as can be verified by the application of theorem 5.1, since *G* satisfies (5.1*a*) and (5.1*b*).

All functions *f*_{ij}(*i*≠*j*), except for {*i*, *j*}∈{{1, 3}, {3, 5}, {2, 4}}, are set equal to *f*^{*}, withFurthermore,The idea is that, for a well-chosen initial condition, *x*_{2} and *x*_{4} remain fixed at −2 and 2, respectively, with agent 3 in between, and agents 1 and 5 going to −∞ and +∞, respectively, and with all agents separated by at least 1 for all *t*≥*t*_{0} (*t*_{0}>0 yet to be defined). Consequently, all interactions are saturated, except for the interaction between agents 1 and 3, 3 and 5, and 2 and 4. In general, the asymptotic relative positions of the agents in a cluster are determined by the internal interactions and the saturation values of external interactions. However, in the present example, the interaction between agents 2 and 3 and agents 3 and 4 are saturated, and the saturation values of the interactions of 1 and 5 on 3 are cancelled. Consequently, the relative position of agent 3 within the cluster {2, 3, 4} depends on the initial condition and on the way the interaction functions *f*_{13} and *f*_{35} approach their saturation values *F*_{13} and *F*_{35}. The dependence on the initial condition has already shown that the result of theorem 5.5 does not hold; we will also choose *f*_{13} and *f*_{35} appropriately to construct a solution to (2.2) for which *x*_{3}(*t*) returns infinitely many times to both −1 and 1, and therefore there is no convergence of *x*_{4}(*t*)−*x*_{3}(*t*) and *x*_{3}(*t*)−*x*_{2}(*t*).

The (asymptotic) relative position of agents 2 and 4 is determined by the choice for the function *f*_{24}, which guarantees that if *x*_{4}(*t*)−*x*_{2}(*t*)=4 for *t*=*t*_{0}, then it remains equal to 4 for *t*>*t*_{0}. Although this choice for *f*_{24} is not necessary and we could have chosen *f*_{24} equal to the function *f*^{*}, it makes sure that (5.1*a*) and (5.1*b*) are satisfied for the cluster structure *G*, showing that this non-generic situation, with distances between agents from the same cluster not converging to a value independent of the initial condition, is not necessarily related to the non-generic situation where (5.2*a*) and (5.2*b*) are satisfied with at least one equality (and therefore (5.1*a*) and (5.1*b*) are not satisfied).

With the interactions between agents 2 and 3 and agents 3 and 4 saturated, the behaviour of agent 3 satisfiesfor all *t*≥*t*_{0}. For *t* sufficiently large, we can assume that *x*_{3}(*t*)−*x*_{1}(*t*) and *x*_{5}(*t*)−*x*_{3}(*t*) are increasing in *t* (as their derivatives approach 1), and therefore *f*_{31}(*x*_{1}(*t*)−*x*_{3}(*t*)) is non-increasing, and *f*_{35}(*x*_{5}(*t*)−*x*_{3}(*t*)) is non-decreasing in *t*. In a first attempt, we can choose *f*_{31} and *f*_{35} as piecewise constant, such that *f*_{31}(*x*_{1}(*t*)−*x*_{3}(*t*)) and *f*_{35}(*x*_{5}(*t*)−*x*_{3}(*t*)) are non-increasing (non-decreasing) and approach their saturation values of −1 (1). When *x*_{3}(*t*)=1, the first term is decreased such that ; when *x*_{3}(*t*)=−1, the second term is increased such that . This results in *x*_{3}(*t*) repeatedly going from −1 to 1 and from 1 to −1. Since *f*_{31} and *f*_{35} approach their saturation values, the derivative will approach zero, and the time needed for *x*_{3}(*t*) to cross the interval [−1,1] will grow unbounded.

In order to obtain a simple expression for *x*_{3}(*t*), we will consider other interaction functions *f*_{31} and *f*_{35} than those proposed in the previous paragraph, with the result based on the same underlying principle. We postulate a deviation for *f*_{31}(*x*_{1}(*t*)−*x*_{3}(*t*)) and *f*_{35}(*x*_{5}(*t*)−*x*_{3}(*t*)) from their saturation values −1 and 1 of the order of *O*(1/*t*) for large *t* values, and we consider the following behaviour for agent 3:(Another option would be to consider *x*_{3}(*t*)=sin(*t*^{α}), with *α*∈(0,1) based on a deviation for *f*_{31}(*x*_{1}(*t*)−*x*_{3}(*t*)) and *f*_{35}(*x*_{5}(*t*)−*x*_{3}(*t*)) from their saturation values of the order of *O*(1/*t*^{1−α}).) As a result, we require thatfor all *t*≥*t*_{0}. For *t* sufficiently large, we can assume that *x*_{3}(*t*)−*x*_{1}(*t*) and *x*_{5}(*t*)−*x*_{3}(*t*) are increasing in *t* (as their derivatives approach 1), and we can put forward the following expressions:(5.3a)and(5.3b)for *t*≥*t*_{0}, where *C*≥0 is to be chosen such that the right-hand sides are non-increasing (non-decreasing) in *t*. As the right-hand side of (5.3*a*) is obviously non-increasing in *t*, we consider the derivative of the right-hand side of (5.3*b*) to *t* and demand thatfor all *t*≥*t*_{0}. We therefore set (since ), and we derive that, for a well-chosen initial condition,for all *t*≥*t*_{0.} To meet the requirements that the arguments of *f*_{31} and *f*_{35} in (5.3*a*) and (5.3*b*) are decreasing (increasing) and that all agents are separated over at least 1, we set . The continuation for the undefined parts of *f*_{31} and *f*_{35} is not important, as long as the functions are non-decreasing and approach their saturation values.

For this solution, the velocities of all agents approach a constant value, but the differences *x*_{3}(*t*)−*x*_{2}(*t*) and *x*_{4}(*t*)−*x*_{3}(*t*) do not converge.

## 6. Outline of the proofs

This section describes the main arguments underpinning the results claimed in §5. For full mathematical details, one is referred to the electronic supplementary material.

### (a) Theorem 5.1

The characteristics of the interaction play a key role in the proof of theorem 5.1. Owing to the anti-symmetry properties of the functions *f*_{ij}, all internal interactions (i.e. interactions between agents in the same cluster) are cancelled when calculating the velocity of the ‘centre of mass’ (weighted with the parameters *γ*_{i}) of a cluster, similar to the cancellation of internal interactions in mechanics. The saturation of the interaction functions implies that the interactions between agents from different clusters can be approximated by their saturation values whenever the agents involved are separated by a sufficiently large distance.

These properties lead to the aforementioned conclusion that, for a solution *x* exhibiting clustering behaviour with respect to *G*=(*G*_{1}, …, *G*_{M}), the asymptotic average cluster velocity depends only on the functions *f*_{ij} through their saturation values *F*_{ij} as follows:for all *k*∈*I*_{M}.

Applying this formula under the assumption of clustering behaviour leads to the approximationfor *t* sufficiently large, and the ordering of the agents and distances growing unbounded with time for agents in different clusters then leads to the conditions (5.2*a*). Since the functions *f*_{ij} are non-decreasing, one similarly derives that(with ⪆ meaning that the right-hand side is an approximate lower bound of the left-hand side), for any two subsets *G*_{k,1} and *G*_{k,2} partitioning *G*_{k}. Since distances between agents from the same cluster remain bounded, this leads to the condition (5.2*b*). This implies the necessity of the inequalities (5.2*a*) and (5.2*b*) for the existence of a solution of (2.2) satisfying clustering behaviour.

For the proof of sufficiency of the conditions (5.1*a*) and (5.1*b*), the main idea is based on the construction of a trapping region *R* for solutions of (2.2): any solution starting in *R* remains in *R* for all later time instances. The trapping region *R* incorporates both a bounded distance between agents of the same cluster and a minimal separation of agents from different clusters. The proof that *R* (defined in the electronic supplementary material) is a trapping region is based on similar expressions as in the previous paragraph in combination with the inequalities (5.1*a*) and (5.1*b*). From the properties of *R*, it then easily follows that any solution *x* starting in *R* satisfies clustering behaviour (with *T*=0 and the clusters equal to the sets *G*_{k}).

Any other solution of (2.2) will exhibit the same clustering behaviour (i.e. identical clusters, possibly a different value for *T*). This follows by observing that we can introduce a modified square distance in the state space between *x* and that is non-increasing in time, due to the monotonicity of the functions *f*_{ij},It follows that remains bounded for all *i* in *I*_{N}, and therefore *x* and exhibit the same clustering behaviour.

### (b) Theorem 5.2

The relationship between the cluster structure and the corresponding intervals for the coupling strength is investigated as follows. The cluster structure *G*=({1}, …, {*N*}) satisfies the conditions (5.2*a*) and (5.2*b*) for *K*=0, and therefore also for values of *K* in some interval [0,*K*_{t}] (with *K*_{t}≥0; if *K*_{t}=0, then (5.2*a*) and (5.2*b*) are satisfied for *K*=0 only—ignoring negative values for *K*). When *K* increases, transitions to a different cluster structure will take place each time one of the inequalities in (5.2*a*) or (5.2*b*) becomes an equality. At the transition value *K*_{t}, a new cluster structure can be constructed, which satisfies (5.2*a*) and (5.2*b*) for values of *K* in some interval (with ). If one of the inequalities in (5.2*a*) becomes an equality at *K*_{t}, then the corresponding clusters *G*_{k} and *G*_{k+1} will merge and form a new cluster. If one of the inequalities in (5.2*b*) becomes an equality at *K*_{t}, then the cluster *G*_{k} will split into two new clusters *G*_{k,1} and *G*_{k,2}, corresponding to the two subsets involved in the equality. (The calculations showing that this new cluster structure satisfies (5.2*a*) and (5.2*b*) for *K* in may be tedious, but they are quite straightforward.) By repeating this procedure, we obtain intervals for the coupling strength *K*, each associated with a particular cluster structure satisfying (5.2*a*) and (5.2*b*) in this interval, and therefore satisfying (5.1*a*) and (5.1*b*) in the interior of this interval. The uniqueness of the cluster structures follows from theorem 5.1.

### (c) Theorem 5.3

To investigate the internal behaviour of a cluster *G*_{k}, *k*∈*I*_{M}, we consider the subsystem that corresponds to the agents in this cluster, with their position relative to the average position of the agents in *G*_{k}, and where the interaction with any agent from another cluster is modelled as the sum of the corresponding saturation value and a time-varying deviation that tends to zero for large *t*, owing to the clustering behaviour. By introducing a Lyapunov function, it is possible to show that (for almost all ) the solution of this subsystem approaches a unique equilibrium point, and therefore that distances between agents from the same cluster will approach a constant value, independent of the initial condition, and that the velocity of an agent converges to the asymptotic average cluster velocity.

## 7. Swarming

As mentioned before, most swarming models focus on the collective behaviour of a species and their cohesion as a single cluster. This may be ensured by imposing periodic or reflective boundary conditions on a bounded region (Czirók & Vicsek 2000), or by including attractive interactions (Gazi & Passino 2003; Mogilner *et al*. 2003). The attraction may decay to zero with increasing mutual distance while still maintaining a coherent cluster (unless it equals zero from some distance on, as in Czirók *et al*. 1996). The size of the population defining a coherent group motivates the short range attraction; indeed, swarms such as schools of fish may consist of several thousands to a million individual animals (Camazine *et al*. 2003). The attraction counteracts the propensity to diverge due to either repulsive terms (Mogilner *et al*. 2003) or drift induced by random walk behaviour (often modelled by inserting a white noise term in the system equations; Mikhailov & Zanette 1999), or both.

In the following paragraphs, we will explain how the model (2.1) may be useful as a model for swarming. For simplicity, we may think of a group of larger land mammals, such as a herd of grazing ungulates.

Since, in the model (2.1), the behaviour of an agent is restricted to one dimension, we interpret its dynamics in applications to swarming as reflecting the behaviour along one direction, thus providing qualitative and not quantitative information. (An extension of the model and some of the results to a multidimensional state space for each agent is possible (De Smet & Aeyels in preparation *a*), but here we focus on the one-dimensional case.)

Other differences between (2.1) and most swarming models described in the literature concern the range of the attraction and the incorporation of the propensity to diverge. In the present model, the range of the attraction is not limited, which for small herds may be a reasonable assumption. Furthermore, the model makes it possible to restrict interactions according to a given network structure, which could be of interest where interaction refers to, for example, kinship.

We now discuss the propensity to diverge, represented by the autonomous components *b*_{i}. They may be interpreted as individual preferences, or as (short- term) average velocities resulting from generalized random walk behaviour (with a continuous range for its increments) due to, for example, external influences. In the latter case, *b*_{i} is time varying on a large time scale, with the cluster structure varying accordingly. Assuming that the clusters develop sufficiently fast, it is possible to observe the cluster structure before another cluster structure emerges. Although, considered on a small time scale, multiple clusters will emerge, they will not persist. The long-term average values for *b*_{i} are small in absolute value (owing to the random walk behaviour, *b*_{i} is of the order of , with *t* denoting time), and (for most time instances) these long-term average values will satisfy the conditions (5.1*a*) and (5.1*b*) for a cluster structure consisting of a single cluster, i.e. *G*=(*I*_{N}). Therefore, the behaviour of the animals will roughly correspond to a single cluster on a larger scale. Animals may drift off to graze further away, but due to the continuous change in cluster structure, they will, on a large scale, form a single cluster; this cluster may not be bounded since some members of the group may continue to wander off in time. (The resulting time-variant model and associated results will be investigated in detail in the future work of De Smet & Aeyels (in preparation *b*).)

### (a) Simulations of the all-to-all coupled model

In our initial simulations, we consider the model from Aeyels & De Smet (2007, 2008), i.e. the system (2.1) with an all-to-all interaction structure, and *A*_{i}=1, for all *i* in *I*_{N}, and with the interaction functions *f*_{ij} all equal to the same function *f*, which is odd, and has a saturation value denoted by *F*. The function *f* is linear in the interval [−2,2] (in m), with slope *F*/2. We set *N*=10, *K*=2 m s^{−1} and *F*=1. For all simulations with *K*=2 m s^{−1}, we use the Euler method with a time step equal to 1 s. Trial showed that smaller time steps were not necessary. The animals start from a single cluster at the position zero: *x*_{i}(0)=0, for all *i*∈*I*_{N}.

In the first simulation, we highlight the emergence of various cluster structures with varying *b*_{i} values by keeping the *b*_{i} values constant during fixed time intervals, allowing us to observe the development of the clusters as predicted by the analytical results. The *b*_{i} values are chosen randomly from a Gaussian distribution with zero mean and a standard deviation of 1 m s^{−1}. They are piecewise constant in time: they are assigned (independently of previous values) every 500 s. Given that the interaction function *f* reaches its saturation value for a distance between interacting animals equal to 2 m, and that the velocities *b*_{i} are of the order of 1 m s^{−1}, we can expect that the time needed for the formation of the clusters is much smaller than the time scale on which the *b*_{i} values vary.

A simulation showing the position of the animals in time is given in figure 2*a*, where the emerging cluster structures can be clearly identified. Given the abrupt changes in position following each discontinuity in *b*_{i} this scenario is unrealistic, except perhaps when modelling the sudden appearance of a predator.

To generate a more natural behaviour in time, in a second simulation (shown in figure 2*b*), the *b*_{i} values are smoother, obtained by sending white noise through a low-pass filter, implemented by the following difference equation in the simulations with a time step of 1 s (with *n* denoting the step number in the Euler integration scheme):(7.1)The random variable *X*_{i}(*n*) is chosen (each second) from a Gaussian distribution with zero mean and a standard deviation of 1 m s^{−1}. The initial value *b*_{i}(0) is drawn randomly from the same Gaussian distribution. The cut-off frequency of the filter is 1/500 Hz, causing a substantial reduction of frequency components with a period smaller than 500 s. The filter dynamics imposed on *b*_{i} result in a Gaussian distribution for *b*_{i}(*n*) (as *b*_{i}(*n*+1) is a linear combination of *b*_{i}(*n*) and *X*_{i}(*n*), which are mutually independent Gaussian random variables). The resulting behaviour, as shown in figure 2*b*, is more natural with respect to the evolution of the cluster structure. Several simulations with the same parameters indicate that a *central cluster* can be distinguished, containing animals with smaller *b*_{i} values (in absolute value). Individual animals or smaller groups, characterized by (temporarily) more extreme *b*_{i} values, occasionally drift off, returning to the main herd later on.

### (b) Simulations with a random network structure

In the next simulation, we consider *N*=50 animals and we abandon the all-to-all interaction structure, but keep *K*=2 m s^{−1} and *A*_{i}=1, for all *i* in *I*_{N}.

First, we choose the network structure for the interactions as a random graph, with each possible edge being selected independently with probability 1/5. The weighting factors are adjusted accordingly, to keep the total attraction experienced by an animal on average equal to the value from the previous simulations, , for all *i* in *I*_{N}. The interaction functions *f*_{ij} are either equal to zero (if there is no link in the interaction network) or equal to the same function *f* as described before, with saturation value *F*=1. The results of a simulation with smooth *b*_{i} values (generated as in the previous simulation) are shown in figure 3.

Again a *central cluster* can be distinguished, with individual animals or smaller groups drifting off and returning to the main herd later on.

### (c) Simulations with a nearest neighbour network structure

As mentioned before, biological models usually consider short-range interactions. In the model (2.1) this may be reflected, without changing the characteristics of the interaction functions, by adopting a nearest neighbour network structure. The network structure is then no longer constant, but is updated each time the order of the *x*_{i}(*t*) values changes. We chose a network structure where each animal interacts with its three nearest neighbours on each side (or less if there are less than three neighbours on one side). We adjusted the weighting factors accordingly, , for all *i* in *I*_{N} (for convenience, no exceptions are made for animals interacting with less than six neighbours). The result is shown in figure 4*a*, and indicates that the interaction is insufficient for the development of clusters.

We therefore increase the coupling strength *K* to 6 m s^{−1}, resulting in figure 4*b*. The simulation time step is decreased to 0.2 s. However, *b*_{i}(*n*) is still updated each second according to (7.1)—with *b*_{i} kept constant between updates—to facilitate a comparison with the previous simulations. In figure 4*b*, one clearly notices the emergence of several clusters on a small spatial scale, which, due to the time variance of the *b*_{i} values and the corresponding cluster structure, constitute a single cluster on a large scale. As opposed to the behaviour in previous simulations (figures 2*b* and 3), this single cluster is not built around a central herd, but consists of several smaller clusters of similar size.

For comparison, the simulation with a random graph is repeated with *K*=6 m s^{−1} (and again a time step of 0.2 s), resulting in figure 5. Here, the increase in coupling strength leads to a single cluster containing all animals, with almost no animals drifting off.

## 8. Compartmental systems

To illustrate the relationship of the proposed clustering model with compartmental systems, we will focus on a system of interconnected water basins with arbitrary connection structure. In Aeyels & De Smet (2008), we considered a similar example, restricted to the all-to-all coupled case. For completeness, we will again introduce and discuss the relationship between the physical system and the clustering model, but now having an arbitrary network structure in mind.

### (a) Interconnected water basins

We consider *N* separate basins connected by pipes, each basin furthermore subject to either a constant external inflow or outflow of water. For laminar flow through a pipe connecting two basins, the fluid velocity is proportional to the pressure difference (Hagen–Poiseuille law). For larger velocities, the flow becomes turbulent and the relationship between velocity and pressure difference is described by the Darcy–Weisbach equation (Plapp 1968),where Δ*p* is the pressure difference; *L* and *D* are the length and diameter of the pipe, respectively; *ρ* is the fluid density; *v* is the mean fluid velocity (i.e. the ratio of the volume flow rate and the cross-sectional area); and *λ* is the friction factor. Although the friction factor depends on the Reynolds number (which is proportional to the fluid velocity), its variation is small for large values of the Reynolds number. We will *approximate* the resulting relationship () by a saturating function, keeping in mind that the flow rate cannot grow unbounded (owing to the finite basin heights and finite resulting pressures). In other words, we assume that the pipes have a maximal throughput, which may depend on the direction of the flow, and is denoted by *F*_{ij} for the pipe connecting basin *i* and *j* in the direction from *j* to *i*, and *F*_{ji} for the other direction. Representing the water height of basin *i* by *x*_{i}, the pressure difference between basins *i* and *j* will be proportional to *x*_{j}−*x*_{i}, and thus the volume flow rate through the connecting pipe can be represented by *f*_{ij}(*x*_{j}−*x*_{i}). The absence of a pipe between basins *i* and *j* corresponds to setting the corresponding saturation values *F*_{ij} and *F*_{ji} equal to zero.

Denoting the inflow for basin *i* by *Q*_{i} and its surface area—which we assume to be water-level independent—by *S*_{i}, one deriveswhich is the model (2.1) with *K*=1, , *γ*_{i}=1 and , for all *i*.

The dynamics will exhibit clustering behaviour as predicted by our analysis, in the sense that the water-level heights of some of the basins will (asymptotically) increase/decrease with the same velocity. The velocity *v*_{k} associated with cluster *G*_{k} can be calculated as

It is clear that the model is valid as long as the basins do not overflow, and all pipes remain below the water level of the basins they connect. We assume that the basins and initial water-level heights are such that these conditions are satisfied during the *transient behaviour*. As explained in the next paragraph, this implies that we can derive the behaviour of the physical system (with respect to overflowing basins or basins running empty) from the cluster structure of the mathematical model and the corresponding long-term velocities, *v*_{k}.

The problem we are interested in is to check whether a network of basins is prone to flooding, i.e. one requires *v*_{k}≤0 for all clusters *G*_{k}. Assume from now on that the total external inflow equals the total external outflow. When, for this case, the dynamical behaviour of the model reveals the existence of more than one cluster there would be at least one cluster *G*_{k} corresponding to flooding, with a positive value for *v*_{k}. Basins will overflow and, after some time, the model will no longer represent the physical situation. However, the model (2.1) continues to be useful for investigating compartmental systems and revealing problematic situations with respect to flooding. Indeed, a simulation of the mathematical model, although only valid in a finite time interval, may reveal, in its long-term behaviour, the existence of multiple clusters; the overflowing basins in the real world system correspond to the clusters with a positive *v*_{k}.

For a solution with two clusters, *G*_{1} and *G*_{2} (with *v*_{1} negative and *v*_{2} positive), the interpretation of the corresponding inequality (5.1*a*) is that the production in the basins belonging to *G*_{2} cannot be transported to the basins belonging to *G*_{1} by their interconnections.

Under the assumption of balanced in- and outflow, a solution with a single cluster *G*_{1}=*I*_{N} has a corresponding velocity equal to . This is a case where no flooding will occur, no basins will run empty and the model will remain valid for all positive time, with each *x*_{i} approaching a constant value.

As an illustration, consider a configuration of *N*=10 basins all having the same surface area *S*_{i}=1, implying that *b*_{i}=*Q*_{i}. (For simplicity we will omit units.) The vector *b* containing the *b*_{i} values is given by(Note that the net inflow in the configuration is zero.) The pipes all have a maximal throughput equal to 2 (in both directions) and are connected in a ring structure, as shown in figure 6*a*. A simulation (figure 7*a*) reveals that this interconnection structure is not able to prevent basins from overflowing (i.e. there are different clusters and at least one of them has a positive velocity). Note that any simulation with arbitrary initial conditions will settle into the same cluster structure, as shown in the proof of theorem 5.1. The objective is to alter the connection structure by adding a minimal number of pipes (of maximal throughput 2) in order to avoid flooding, i.e. in order to obtain a single cluster at zero velocity.

Figure 7*a* shows that there are six different clusters: *G*_{1}={1}; *G*_{2}={4, 5}; *G*_{3}={8, 9, 10}; *G*_{4}={6}; *G*_{5}={2, 3}; and *G*_{6}={7}. Adding an extra pipe will not affect the velocity of the cluster with largest (smallest) *v*_{k} unless the pipe is connected to one of the basins in this cluster. Therefore, if one extra connection would be sufficient it would have to connect basins 1 and 7. A simulation with this extra connection results in three clusters: ; ; and (figure 7*b*), implying that we still need (at least) one extra connection between a basin from and . Simulations show that any extra connection between either 4 or 5 and 2 or 3 leads to one cluster at zero velocity, solving the problem. A possible solution is shown in figure 6*b*.

This example shows that the model (2.1) is useful for the analysis of compartmental systems and in suggesting solutions to associated problems.

### (b) Other examples

Other applications of compartmental systems may require specific adaptations of the model and may suggest extensions.

Consider a system of lakes interconnected by rivers or channels, for which the lakes, rivers and channels are not allowed to overflow. The volume flow rate of rivers emanating from

*artificial*lakes can be controlled and, accordingly, one can prevent them from flooding. The main objective is then to avoid the reservoirs from overflowing. Rivers emanating from*natural*lakes have no control mechanism to restrict the water supply to their maximal throughput, imposing a second objective of preventing the*rivers*from overflowing. In both cases, the volume flow rate of a river flowing from one lake to another will not be a mere function of the water-level difference between the lakes (for natural lakes, the volume flow rate mainly depends on the absolute water-level height of the lake), and the model (2.1) will not be able to describe the behaviour of the water-level heights. However, the inequalities (5.2*b*) express that the capacity of the channels and the rivers has to be sufficiently large for transporting the excess water out of (the lakes in) regions with a positive net supply rate. Therefore, these inequalities are*necessary*to prevent flooding. They can be checked directly or by simulation of the differential equations (2.1). To obtain a set of*sufficient*conditions for flood prevention, additional conditions (e.g. to include directionality of connections or expressing that rivers should not overflow) need to be included.A major issue regarding road traffic is to avoid congestions. Different regions can be considered as compartments interconnected by highways, which have limited transport capacities. Again, the expression for the flow rate needs to be modified, but the main property—the saturation—still holds. Regarding regular traffic, the model will need to incorporate the fact that cars are not interchangeable since each driver has its own destination. This aspect can be neglected in special cases (e.g. in the holiday season), where there is a common destination for the majority of the drivers and then the remark on necessity of the inequalities in the previous item still applies.

## 9. Opinion formation

We represent opinions on a particular matter by real numbers, with zero corresponding to a neutral position. We consider *N* individuals taking part in a meeting; each individual has his own opinion on the issue on the agenda, which may evolve in time due to discussion with the other members. Since opinions cannot grow unbounded, *x*_{i} in (2.1) is not an appropriate quantity to represent an opinion. Instead, we will take the derivatives as a measure of someone's opinion. The equations for *y*_{i} can be written as (assuming *x*_{i}(0)=0, ∀*i*∈*I*_{N}, without loss of generality regarding the long-term behaviour)(9.1)∀*i*∈*I*_{N}, where we have redefined the sensitivity factors *A*_{i} to explicitly include a normalization of the interaction, such that each agent's opinion deviates at most *KA*_{i} from its *a priori* value *b*_{i} (corresponding to ‘no discussion’). With *y*_{i}(*t*) representing the opinion of agent *i* at time *t*, the absolute value of the integral may reflect the level of disagreement accumulated over time, or the amount of discussions taking place between agents *i* and *j*, proportional with time and difference in opinion. We assume that all *F*_{ij} are non-negative, and that *f*_{ij}(0)=0 for all *i*, *j*∈*I*_{N}; interactions are attractive and take place only if there is discussion.

In general, everyone starts with his own opinion *b*_{i}, while, with time and through interaction, different groups are formed, each group is characterized by a final opinion *v*_{i} obtained through discussion. The pressure to reach a decision or the tendency to adapt one's opinion by paying attention to each other's arguments is reflected by the coupling strength *K*.

In figure 8*a*, we show the evolution of the opinions *v*_{i} eventually reached as a function of *K*. (The other parameters are left unchanged.) The *v*_{i} values were calculated by means of an algorithm based on the inequalities (5.1*a*) and (5.1*b*), not by a simulation of the integral equation (9.1). We considered 100 agents with *b*_{i} chosen from a Gaussian distribution with zero mean and a standard deviation of 1. The parameters *A*_{i}, *γ*_{i} and *F*_{ij} (*i*≠*j*) were all taken to be equal to 1. Notice a steady convergence to complete agreement as a function of *K*. In figure 8*b*, the time evolution of the opinions *y*_{i} for *K*=1.5, as obtained by numerical integration of the mathematical model, is shown.

In a second simulation (figure 8*c*), we kept the same parameters *b*_{i}, but the values for *A*_{i} and *γ*_{i} were changed to account for the fact that people with extreme opinions are reluctant to change their point of view (smaller *A*_{i}), while making more efforts to persuade other people (larger *γ*_{i}). Also, *F*_{ij} decreases with increasing values of |*b*_{j}−*b*_{i}|, reflecting the idea that people tend to pay more attention to people with similar opinions,(9.2)(9.3)and(9.4)In figure 8*d*, we show the time evolution of *y*_{i} for *K*=15, again obtained by numerical integration. (For the numerical integration in figure 8*b*,*d* the Euler method was used with a time step of 0.03/*K*.)

While, in the first case, it seems possible to take a decision by unanimous consent, in the second case—which is more realistic—it is far more favourable to let a majority vote decide, as one notices a deadlock of extreme opinions for a *K* of approximately 15. Total consensus can only be reached under much higher pressure compared with the pressure needed for reaching a decision by a majority vote, and might require unreasonable concessions from all parties involved.

As an important distinction with other existing models (for an overview, see Hegselmann & Krause 2002), we want to emphasize that the model (9.1) allows the coexistence of several groups, each characterized by its own group opinion—as opposed to models focusing on total consensus or the coexistence of only two opinions (such as in Sznajd-Weron & Sznajd 2000)—while still allowing analytical exploration, as opposed to models for which the results rely on simulations (as in Deffuant *et al*. 2002) or statistical methods (e.g. Bordogna & Albano 2007).

## Acknowledgments

We thank the anonymous referees for their constructive comments.

This paper presents research results of the Belgian Programme on Interuniversity Attraction Poles, initiated by the Belgian Federal Science Policy Office. The scientific responsibility rests with its authors.

During this research, F.D.S. was supported by a PhD fellowship of the Research Foundation—Flanders (FWO).

## Footnotes

- Received June 19, 2008.
- Accepted October 23, 2008.

- © 2008 The Royal Society