Tiziano Squartini

Diego Garlaschelli

Triadic motifs and dyadic self-organization in the World Trade Network

In self-organizing networks, topology and dynamics coevolve in a continuous feedback, without exogenous driving. The World Trade Network (WTN) is one of the few empirically well documented examples of self-organizing networks: its topology depends on the GDP of world countries, which in turn depends on the structure of trade. Therefore, understanding the WTN topological properties deviating from randomness provides direct empirical information about the structural effects of self-organization. Here, using an analytical pattern-detection method we have recently proposed, we study the occurrence of triadic 'motifs' (three-vertices subgraphs) in the WTN between 1950 and 2000. We find that motifs are not explained by only the in-and out-degree sequences, but they are completely explained if also the numbers of reciprocal edges are taken into account. This implies that the self-organization process underlying the evolution of the WTN is almost completely encoded into the dyadic structure, which strongly depends on reciprocity.

Introduction

The global economy is a prototypic example of complex self-organizing system, whose collective properties emerge spontaneously through many local interactions. In particular, international trade between countries defines a complex network which arises as the combination of many independent choices of firms. It was shown that the topology of the World Trade Network (WTN) strongly depends on the Gross Domestic Product (GDP) of world countries [START_REF] Garlaschelli | Fitness-dependent topological properties of the World Trade Web[END_REF]. On the other hand, the GDP depends on international trade by definition [START_REF] Garlaschelli | Interplay between topology and dynamics in the World Trade Web[END_REF], which implies that the WTN is a remarkably well documented example of adaptive network, where dynamics and topology coevolve in a continuous feedback. In general, understanding self-organizing networks is a major challenge for science, as only few models of such networks are analytically solvable [START_REF] Garlaschelli | Self-organized network evolution coupled to extremal dynamics[END_REF]. However, in the particular case of the WTN, the binary topology of the network is found to be extremely well reproduced by a null model which incorporates the degree sequence [START_REF] Squartini | Randomizing world trade. I. A binary network analysis[END_REF]. These results, which have been obtained using a fast network randomization method that we have recently proposed [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF], make the WTN particularly interesting. In this paper, after briefly reviewing our randomization method, we apply it to study the occurrence of triadic 'motifs', i.e. directed patterns involving three vertices (see Fig. 1). We show that, unlike other properties which have been studied elsewhere [START_REF] Squartini | Randomizing world trade. I. A binary network analysis[END_REF], the occurrence of motifs is not explained by only the in-and out-degrees of vertices. However, if also the numbers of reciprocal links of each vertex (the reciprocal degree sequence) are taken into account, the occurrences of triadic motifs are almost completely reproduced. This implies that, if local information is enhanced in order to take into account the reciprocity structure, motifs display no significant deviations from random expectations. Therefore the (in principle complicated) self-organization process underlying the evolution of the WTN turns out to be relatively simply encoded into the local dyadic structure, which separately specifies the number of reciprocated and non-reciprocated links per vertex. Thus the dyadic structure appears to carry a large amount of information about the system.

Searching for Non-Random Patterns in the WTN

In this section we briefly summarize our recently proposed randomization method and how it can be used to detect patterns when local constraints are considered.

Network Pattern Detection: the Randomization Method

Our method, which is based on the maximum-likelihood estimation of maximumentropy models of graphs, introduces a family of null models of a real network and uses it to detect topological patterns analytically [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF]. Defining a null model means setting up a method to assign probabilites [START_REF] Maslov | Detection of topological patterns in complex networks: correlation profile of the Internet[END_REF][START_REF] Holland | Sociological Methodology[END_REF][START_REF] Park | Statistical mechanics of networks[END_REF]. In our approach, a real network G * with N vertices is given (either a binary or a weighted graph, and either directed or undirected, whose generic entry is g ij ) and a way to generate a family G of randomized variants of G * is provided, by assigning each graph G ∈ G a probability P (G). In the method, the probabilities P (G) are such that a maximally random ensemble of networks is generated, under the constraint that, on average, a set {C a } of desired topological properties is set equal to the values {C a (G * )} observed in the real network G * . This is achieved as the result of a constrained Shannon-Gibbs entropy maximization [START_REF] Park | Statistical mechanics of networks[END_REF] S ≡ -

∑ G P (G) ln P (G) (1)
and the imposed constraints are the normalization of the probability and the average values of a number of chosen properties, {C a }:

1 = ∑ G P (G), ⟨C a ⟩ ≡ ∑ G C a (G)P (G) . ( 2
)
This optimization leads to exponential probability coefficients [START_REF] Park | Statistical mechanics of networks[END_REF]

P (G) = e -H(G,θ) Z (3)
where the linear combination of the contraints H(G, θ) ≡ ∑ a θ a C a (G) is called graph hamiltonian (the coefficients {θ a } are free parameters, acting as Lagrange multipliers controlling the expected values {⟨C a ⟩}) and the denominator Z ≡ ∑ G e -H(G,θ) is called partition function. The next step is the maximization of the probability P (G * ) to obtain the observed graph G * , i.e. the real-world network to randomize [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF]. This step fixes the values of the Lagrange multipliers as they are found by maximizing the log-likelihood

L(θ) ≡ ln P (G * |θ) = -H(G * , θ) -ln Z(θ) (4)
to obtain the real network G * . It can be easily verified [START_REF] Garlaschelli | Maximum likelihood: extracting unbiased information from complex networks[END_REF] that this is achieved by the parameter values θ * satisfying

⟨C a ⟩ * = ∑ G C a (G)P (G|θ * ) = C a (G * ) ∀a (5)
that is, that the ensemble average of each constraint, ⟨C a ⟩, equals the observed value on the real network, C a (G * ). Once the numerical values of the Lagrange multipliers are found, they can be used to find the ensemble average ⟨X⟩ * of any topological property X of interest [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF]:

⟨X⟩ * = ∑ G X(G)P (G|θ * ) . ( 6
)
The exact computation of the expected values can be very difficult. For this reason it is often necessary to rest on the linear approximation method [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF]. However, in the present study we will consider particular topological properties X (i.e. motif counts, see below) whose expected value can be evaluated exactly. Our method also allows to obtain the variance of X by applying the usual definition:

σ 2 [X] = ⟨[X(G) -⟨X⟩)] 2 ⟩ = ∑ i,j ∑ t,s σ[g ij , g ts ] ( ∂X ∂g ij ∂X ∂g ts ) G=⟨G⟩ (7)
where σ[g ij , g ts ] is the covariance of the adjacency matrix elements g ij and g ts . This formula can be greatly simplified by considering probabilities that are factorizable in terms of dyadic probabilities, as follows

P (G|θ) = ∏ i<j D ij (g ij , g ji |θ) ( 8
)
where the product runs over all the dyads, that is the unordered pairs of vertices (i, j) (with i < j), and D ij (g, g ′ |θ) is the joint probability that g ij = g and g ji = g ′ . Finally, the variance of X evaluated in θ * becomes

(σ * [X]) 2 = ∑ i,j [ ( σ * [g ij ] ∂X ∂g ij ) 2 G=⟨G⟩ * + σ * [g ij , g ji ] ( ∂X ∂g ij ∂X ∂g ji ) G=⟨G⟩ * ] . ( 9
)
The joint knowledge of ⟨X⟩ * and σ * [X] allows to detect deviations from randomness in the observed topology. In particular, as we show later, it is possible to calculate by how many standard deviations the observed value X * differs from the expected value ⟨X⟩ * . Quantities which are consistent with their expected value are explained by the enforced constraints {C a }. On the other hand, significantly deviating properies cannot be traced back to the constraints and therefore signal the incompleteness of the information encoded in the constraints. Other approaches achieve this result by explicitly generating many randomized variants of the real network, measuring X on each such variant, and finally computing the sample average and standard deviation of X [START_REF] Maslov | Detection of topological patterns in complex networks: correlation profile of the Internet[END_REF]. This is extremely time consuming, especially for complicated topological properties. By contrast, our method is entirely analytical. It yields any expected quantity ⟨X⟩ * in a time as short as that required in order to measure X * on the single network G * [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF].

The Role of Local Constraints

If the network is a binary graph (i.e. if each graph G in the ensemble is uniquely specified by its adjacency matrix A), then the simplest (i.e. local) choice of the constraints is the degree sequence, i.e. the vector of degrees (numbers of incident links) of all vertices. For directed networks, which are our interest here, there are actually two degree sequences: the observed in-degree sequence k in (A * ) (with k in i = ∑ j̸ =i a ji ) and the observed out-degree sequence k out (A * ) (with

k out i = ∑ j̸ =i a ij
). This null model, which is known as the directed configuration model (DCM), can be completely dealt with analytically using our method (see Appendix). When applied to the WTN, the DCM shows that many topological properties (such as the degree-degree correlations and the directed clustering coefficients) are in complete accordance with the expectations [START_REF] Squartini | Randomizing world trade. I. A binary network analysis[END_REF]. This shows that the degree sequences k in and k out are extremely informative, as their (partial) knowledge allows to reconstruct many aspects of the (complete) topology. On the other hand, it was also shown that the reciprocity of the WTN is highly non-trivial [START_REF] Garlaschelli | Patterns of link reciprocity in directed networks[END_REF]. This means that the occurrence of reciprocal links is much higher than expected under any model which, as the DCM, treats two reciprocal links (e.g. i → j and j → i) as statistically independent. A direct consequence is that the reciprocity, as well as any higher-order directional pattern, should not be reproduced by the DCM. These seemingly conflicting results can only be reconciled if, for some reason, the topological properties that have been studied under the DCM [START_REF] Squartini | Randomizing world trade. I. A binary network analysis[END_REF] mask the effects of reciprocity. In particular, the directed clustering coefficients, which are based on ratios of realized triangles over the maximum number for each vertex, may show no overall deviation from the DCM, even if the numerator and denominator separately deviate from it. In what follows, we investigate this possibility by considering all the observed subgraphs of three vertices (which include both open and closed triangles) separately. Also, we will use an additional null model which also takes the number of reciprocal links of each vertex into account. This second null model is the reciprocal configuration model (RCM) [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF][START_REF] Garlaschelli | Patterns of link reciprocity in directed networks[END_REF][START_REF] Garlaschelli | Multispecies grand-canonical models for networks with reciprocity[END_REF]. The local constraints defining it are the three, observed directed-degree sequences k → (A * ), with

k → i ≡ ∑ j̸ =i a → ij (non-reciprocated out degree), k ← (A * ), with k ← i ≡ ∑ j̸ =i a ← ij (non-reciprocated in-degree) and k ↔ (A * )}, with k ↔ i ≡ ∑ j̸ =i a ↔ ij (reciprocated degree
) to be imposed across the ensemble of networks having the same number of vertices of the observed configuration and, on average, the above-mentioned directed-degree sequences. In the Appendix we describe both the DCM and the RCM in more detail, and derive their expectations explicitly.

Triadic Motifs in the WTN

In the following analyses, we use yearly bilateral data on exports and imports from the Gleditsch Database4 to analyse the six years 1950,1960,1970,1980,1990,2000. This database contains aggregated trade data between countries, i.e. data as they result by summing the single commodity-specific trade exchanges. So we end up with six different, real, asymmetric matrices with entries agg ij (y) (y = 1950, 1960 . . . 2000). These adjacency matrix elements are the fundamental data allowing us to obtain all the possible representations of the WTN: to build the binary, directed representation we are interested in here, we restrict ourselves to consider two different vertices as linked, whenever the corresponding element m agg ij (y) is strictly positive. This implies that the adjacency matrix of the binary, directed representation of the WTN in year y is simply obtained by applying the Heaviside step function to the database entries, i.e. a ij (y) = Θ[m agg ij (y)]. Triadic motifs, i.e. all the possibile directed patterns connecting three vertices, are the natural generalizations of directed clustering coefficients [START_REF] Fagiolo | Clustering in complex directed networks[END_REF] and the starting point for the understanding of a complex network's self-organization in communities. Thirteen, non-isomorphic, triadic directed patterns can be indentified and classified [START_REF] Caldarelli | Scale-free Networks[END_REF]. Given a real, binary, directed matrix A * , the motifs occurrences N m can be written in at least two different ways (see Table 1). The first one prescribes to define them in terms of the adjacency matrix entries, {a ij }. The second one allows to compactly express the coefficients N m by introducing the following dyadic variables

a → ij ≡ a ij (1 -a ji ), a ← ij ≡ a ji (1 -a ij ), a ↔ ij ≡ a ij a ji , a ij ≡ (1 -a ij )(1 -a ji ) (10)
thus making the role of reciprocity explicit. However, the number of occurrences of the particular motif m (where m ranges from 1 to 13) as measured on the observed network A * is uninformative unless a comparison with a properly defined null model is made (see Appendix). This implies that the occurrence of a motif should be compared with its expected value ⟨N m ⟩ * , as computed under the chosen null model. This can be compactly achieved by introducing the so-called z-score Table 1. Classification (after [START_REF] Caldarelli | Scale-free Networks[END_REF]) and definitions of the triadic motifs Motif m Nm:

z[N m ] ≡ N m (A * ) -⟨N m ⟩ * σ * [N m ] (11)
1 st definition Nm: 2 nd definition 1 ∑ i̸ =j̸ =k (1 -aij)ajia jk (1 -a kj )(1 -a ik )(1 -a ki ) ∑ i̸ =j̸ =k a ← ij a → jk a ik 2 ∑ i̸ =j̸ =k aij(1 -aji)a jk (1 -a kj )(1 -a ik )(1 -a ki ) ∑ i̸ =j̸ =k a → ij a → jk a ik 3 ∑ i̸ =j̸ =k aijajia jk (1 -a kj )(1 -a ik )(1 -a ki ) ∑ i̸ =j̸ =k a ↔ ij a → jk a ik 4 ∑ i̸ =j̸ =k (1 -aij)(1 -aji)a jk (1 -a kj )a ik (1 -a ki ) ∑ i̸ =j̸ =k a ij a → jk a → ik 5 ∑ i̸ =j̸ =k (1 -aij)ajia jk (1 -a kj )a ik (1 -a ki ) ∑ i̸ =j̸ =k a ← ij a → jk a → ik 6 ∑ i̸ =j̸ =k aijajia jk (1 -a kj )a ik (1 -a ki ) ∑ i̸ =j̸ =k a ↔ ij a → jk a → ik 7 ∑ i̸ =j̸ =k aijaji(1 -a jk )a kj (1 -a ik )(1 -a ki ) ∑ i̸ =j̸ =k a ↔ ij a ← jk a ik 8 ∑ i̸ =j̸ =k aijajia jk a kj (1 -a ik )(1 -a ki ) ∑ i̸ =j̸ =k a ↔ ij a ↔ jk a ik 9 ∑ i̸ =j̸ =k (1 -aij)aji(1 -a jk )a kj a ik (1 -a ki ) ∑ i̸ =j̸ =k a ← ij a ← jk a → ik 10 ∑ i̸ =j̸ =k (1 -aij)ajia jk a kj a ik (1 -a ki ) ∑ i̸ =j̸ =k a ← ij a ↔ jk a → ik 11 ∑ i̸ =j̸ =k aij(1 -aji)a jk a kj a ik (1 -a ki ) ∑ i̸ =j̸ =k a → ij a ↔ jk a → ik 12 ∑ i̸ =j̸ =k aijajia jk a kj a ik (1 -a ki ) ∑ i̸ =j̸ =k a ↔ ij a ↔ jk a → ik 13 ∑ i̸ =j̸ =k aijajia jk a kj a ik a ki ∑ i̸ =j̸ =k a ↔ ij a ↔ jk a ↔

Results and Discussion

Fig. 2 and Fig. 3 show the z-scores for all the 13, triadic, binary, directed motifs, computed for the six different snapshots of the World Trade Network corresponding to the decades 1950,1960,1970,1980,1990,2000, for both the DCM and the RCM. We also show the six lines z = ±1, z = ±2 and z = ±3 to highlight the region within 3 sigmas from the expectation value. The analysis reveals a dramatic difference between the predictions of the two null models. The presence of intrinsically directed trading relationships implies that reciprocity is a fundamental quantity, shaping the network of exchanges among world-countries. The reciprocity of the WTN is known to be very high [START_REF] Garlaschelli | Patterns of link reciprocity in directed networks[END_REF] and this has strong effects on its motif structure. This implies that, in order to reproduce the topology of the network, it is essential to reproduce its dyadic structure, which in turn strongly depends on reciprocity. This also implies that the RCM should fluctuate less than the DCM, because of this huge amount of additional information. Indeed, this is confirmed by our analysis.

In fact, we find that, unlike other topological properties [START_REF] Squartini | Randomizing world trade. I. A binary network analysis[END_REF], triadic motifs are systematically not reproduced if only the in-degree and out-degree sequences are taken into account. In particular, the observed motifs counts deviate by up to 100 standard deviations from the expected ones. By contrast, if also the reciprocity (which separately specifies the number of reciprocated links) is introduced in the null model, triadic motifs are almost always consistent with expectations within one standard deviation, as for m = 5, 10, 11, 12 or at most two standard deviations, as m = 1, 4. Moreover, the significance profiles are almost completely inverted: some of the motifs being over(under)-represented without the reciprocity constraint, as m = 8 (m = 2), become under(over)-represented with the reciprocity constraint.

Conclusions

The WTN is a particularly interesting network which is known to be driven by a self-organization process involving the global economy. Empirically, it is a very well documented network which allows to test predictions of null models about the key ingredients shaping its topology. In this paper, we aimed at isolating the key properties of the WTN topology where the effects of self-organizations can be clearly detected as deviations from randomness. While the DCM, which takes into account the number of incoming and outgoing connections of all vertices, reproduces many other topological properties, it cannot replicate the observed triadic motifs. On the other hand, we found that the RCM, which also takes into account the numbers of reciprocated links, can replicate almost perfectly the triadic structure. We therefore found that the process underlying the evolu-tion of the WTN is mainly encoded into the dyadic structure, which carries a large amount of information about the system: so, the upgrade to the RCM is necessary and the possibility to treat the RCM analytically using our method is therefore an important step forward. The result that dyadic properties almost completely explain triadic ones suggests that in the WTN also higher-order properties (e.g. subgraphs of four or more vertices, and in general even the existence of denser communities of vertices) will be explained mostly by dyadic properties. This raises the important question whether the available community detection methods are successful in identifying communities which are not explained by local constraints. In particular, our results suggest that modularity-based community detection methods will detect communities in the WTN if a null model without reciprocity (DCM) is used, while they should find weak or no community structure if a null model including reciprocity (RCM) is used. In both cases, as we have already claimed [START_REF] Squartini | Analytical maximum-likelihood method to detect patterns in real networks[END_REF], the available expressions for the modularity (which make the strong assumption of sparse networks) should be revised using the correct expectation values provided by our method, since the WTN is a very dense network. We will address these issues in future work.

Appendix

The Directed Configuration Model

Given a real, binary, directed network A * with out-degree and in-degree sequences k out (A * ) and k in (A * ), the method described in section 2.1 can be specified in the following way.

The DCM Hamiltonian. The Hamiltonian implementing the DCM is

H(A, α, β) = ∑ i [α i k out i (A) + β i k in i (A)] = ∑ i̸ =j (α i + β j )a ij ; (12)
the partition function can be calculated as in [START_REF] Park | Statistical mechanics of networks[END_REF], Z(α, β) = ∑ A e -H(A,α,β) = ∏ i̸ =j (1 + e -αi-βj ), and the graph probability is

P (A|α, β) ≡ ∏ i<j D ij (a ij , a ji |α i , α j , β i , β j ) = ∏ i̸ =j p aij ij (1 -p ij ) (1-aij ) (13)
and, by setting x i ≡ e -αi and y i ≡ e -βi ,

p ij = e -αi-βj 1 + e -αi-βj ≡ x i y j 1 + x i y j . ( 14
)
The Log-Likelihood Equations. The log-likelihood function to maximize is

L(x, y) = ∑ i [ k out i (A * ) ln x i + k in i (A * ) ln y i ] - ∑ i̸ =j ln(1 + x i y j ) (15)
and the values x * , y * corresponding to the point of maximum can be found by solving the following system

   k in i (A * ) = ∑ j̸ =i x * i y * j 1+x * i y * j ∀i k out i (A * ) = ∑ j̸ =i x * j y * i 1+x * j y * i ∀i ( 16
)
Expectation Values and Variances. The expectation value and variance of a ij are, respectively,

⟨a ij ⟩ * = p * ij and (σ * [a ij ]) 2 = ⟨a ij ⟩ * (1 -⟨a ij ⟩ * ) = p * ij (1 -p * ij )
. Distinct pairs of vertices are independent random variables; since the first definition of the motifs occurrences only involves products of the adjacency matrix entries, their expectation values can be easily computed, as shown in Table 2. The variance of N m becomes

(σ * [N m ]) 2 = ∑ i,j [ ( σ * [a ij ] ∂N m ∂a ij ) 2 A=⟨A⟩ * ] (17) considering that σ * [a ij , a ji ] = ⟨a ij a ji ⟩ * -⟨a ij ⟩ * ⟨a ji ⟩ * = 0.

The Reciprocal Configuration Model

Given a real binary directed network A * with directed-degree sequences k → (A * ), k ← (A * ) and k ↔ (A * ) where

k → i (A * ) ≡ ∑ j̸ =i a * ij (1 -a * ji ), k ← i (A * ) ≡ ∑ j̸ =i a * ji (1 -a * ij ), k ↔ i (A * ) ≡ ∑ j̸ =i a * ij a * ji ( 18
) the randomization procedure can be specified in the following way.

The RCM Hamiltonian. The Hamiltonian implementing the RCM is

H(A, α, β, γ) = ∑ i [α i k → i (A) + β i k ← i (A) + γ i k ↔ i (A)] , (19)
the partition function can be calculated as in [START_REF] Garlaschelli | Multispecies grand-canonical models for networks with reciprocity[END_REF], Z(α, β, γ) = ∏ i<j (1 + e -αi-βj + e -αj -βi + e -γi-γj ), and the graph probability is

P (A|α, β, γ) = ∏ i<j D ij (a ij , a ji |α i , α j , β i , β j , γ i , γ j ) = ∏ i<j (p → ij ) a → ij (p ← ij ) a ← ij (p ↔ ij ) a ↔ ij (p ij ) a ij (20)
and, by setting x i ≡ e -αi , y i ≡ e -βi and z i ≡ e -γi [START_REF] Garlaschelli | Multispecies grand-canonical models for networks with reciprocity[END_REF],

p → ij ≡ x i y j 1 + x i y j + x j y i + z i z j , p ← ij ≡ x j y i 1 + x i y j + x j y i + z i z j (21) p ↔ ij ≡ z i z j 1 + x i y j + x j y i + z i z j , p ij ≡ 1 1 + x i y j + x j y i + z i z j (22)
The Log-Likelihood Equations. The log-likelihood function to maximize is

L(x, y, z) = ∑ i [k → i (A * ) ln x i + k ← i (A * ) ln y i + k ↔ i (A * ) ln z i ] - ∑ i<j ln(1 + x i y j + x j y i + z i z j ) (23)
and the values x * , y * , z * corresponding to the point of maximum can be found by solving the following system

         k → i (A * ) = ∑ j̸ =i x * i y * j 1+x * i y * j +x * j y * i +z * i z * j ∀i k ← i (A * ) = ∑ j̸ =i x * j y * i 1+x * i y * j +x * j y * i +z * i z * j ∀i k ↔ i (A * ) = ∑ j̸ =i z * i z * j 1+x * i y * j +x * j y * i +z * i z * j ∀i ( 24
)
Expectation Values and Variances. The expectation value and variance of a → ij (and equivalently for the other dyadic variables) are, respectively, ⟨a

→ ij ⟩ * = (p → ij ) * and (σ * [a → ij ]) 2 = ⟨a → ij ⟩ * (1 -⟨a → ij ⟩ * ) = (p → ij ) * (1 -(p → ij ) * ).
Considering that two distinct dyads can be treated as independent random variables and that the second definition of the motifs occurrences only involves products of dyads, their expectation values can be easily computed, as shown in Table 2

Fig. 1 .

 1 Fig. 1. The triadic, binary, directed motifs

ik

  measuring by how many standard deviations, σ * [N m ], the observed and the expected occurrences of motif m differ. Large absolute values of z[N m ] indicate motifs that are either over-or under-represented under the particular null model considered and therefore not explained by the constraints defining it, as shown in Fig.2and Fig.3and discussed in the next section.

Fig. 2 .

 2 Fig. 2. z-scores of the 13 triadic, binary, directed motifs for the six decades of the WTN, under the DCM ( ) and under the RCM (•). The dashed, red lines represent the values z = ±3, the dotted, purple lines the values z = ±2 and the dot-dashed, pink lines the values z = ±1.

Fig. 3 .

 3 Fig. 3. z-scores of the 13 triadic, binary, directed motifs for the six decades of the WTN, under the DCM ( ) and under the RCM (•): zoom of Fig. 2. The dashed, red lines represent the values z = ±3, the dotted, purple lines the values z = ±2 and the dot-dashed, pink lines the values z = ±1.

Table 2 .

 2 . The variance of N m becomes (σ * [N m ]) 2 = Expectation values of the triadic motifs Motif m ⟨Nm⟩DCM ⟨Nm⟩RCM1 ∑ i̸ =j̸ =k (1 -pij)pjip jk (1 -p kj )(1 -p ik )(1 -p ki ) ∑ i̸ =j̸ =k p ← ij p → jk p ik 2 ∑ i̸ =j̸ =k pij(1 -pji)p jk (1 -p kj )(1 -p ik )(1 -p ki ) ∑ i̸ =j̸ =k p → ij p → jk p ik 3 ∑ i̸ =j̸ =k pijpjip jk (1 -p kj )(1 -p ik )(1 -p ki ) ∑ i̸ =j̸ =k p ↔ ij p → jk p ik 4 ∑ i̸ =j̸ =k (1 -pij)(1 -pji)p jk (1 -p kj )p ik (1 -p ki ) ∑ i̸ =j̸ =k p ij p → jk p → ik 5 ∑ i̸ =j̸ =k (1 -pij)pjip jk (1 -p kj )p ik (1 -p ki ) i̸ =j̸ =k pijpjip jk (1 -p kj )p ik (1 -p ki ) -p jk )p kj (1 -p ik )(1 -p ki ) ∑ i̸ =j̸ =k p ↔ ij p ← jk p ik 8 ∑ i̸ =j̸ =k pijpjip jk p kj (1 -p ik )(1 -p ki ) ∑ i̸ =j̸ =k p ↔ ij p ↔ jk p ik 9 ∑ i̸ =j̸ =k (1 -pij)pji(1 -p jk )p kj p ik (1 -p ki ) -pij)pjip jk p kj p ik (1 -p ki ) -pji)p jk p kj p ik (1 -p ki ) ∑ i̸ =j̸ =k p → ij p ↔ jk p → ik 12 ∑ i̸ =j̸ =k pijpjip jk p kj p ik (1 -p ki ) ∑ i̸ =j̸ =k p ↔ ij p ↔jk p → ik 13 ∑ i̸ =j̸ =k pijpjip jk p kj p ik p ki ∑ i̸ =j̸ =k p ↔ ij p ↔ jk p ↔

	∑ i,j	[ ( σ * [a ij ]	∂N m ∂a ij	) 2 A=⟨A⟩ (	∂a ij ∂N m	∂a ji ∂N m	A=⟨A⟩ * )	]
							(25)
	where now (σ

* + σ * [a ij , a ji ] * [a ij ]) 2 = ⟨a ↔ ij + a → ij ⟩ * (1 -⟨a ↔ ij + a → ij ⟩ * ) and σ * [a ij , a ji ] = ⟨a ↔ ij ⟩ * -⟨a ↔ ij + a → ij ⟩ * ⟨a ↔ ji + a → ji ⟩ * .

Gleditsch, K.S.: Expanded trade and GDP data. Jour. Confl. Res. 46 (2002) 712-724.

Acknowledgements

This work was supported by the Dutch Econophysics Foundation (Stichting Econophysics, Leiden, The Netherlands) with funds from beneficiaries of Duyfken Trading Knowledge BV, Amsterdam, The Netherlands.


Laurent Paquereau
email: laurent@item.ntnu.no

Bjarne E Helvik
email: bjarne@q2s.ntnu.no

Ant-Based Systems for Wireless Networks: Retrospect and Prospects

Since they were first introduced as powerful stochastic optimization systems, ant-based systems have been applied to a wide range of problems. One of the most successful applications has been routing in dynamic wired telecommunication networks. Following this success, similar approaches have been applied to routing in multi-hop wireless networks. The objective has been to achieve self-organizing path management in these networks. This paper looks back on 10 years of research and presents reflections on the challenges, the evolution, the contributions and the future perspectives in this field.

Introduction

Ant-based systems are systems inspired by the foraging behavior of ants in nature or, more generally, systems that can be described as such. Ant-based systems are swarm intelligence systems. Each individual ant is not capable of finding the shortest path between its nest and a food source. It is by their collective behavior that ants in a colony are able to converge to the shortest path; a key enabler to this emergent behavior being the indirect communication between ants through the modification of their environment (stigmergy).

Ant-based systems were first developed as powerful stochastic optimization systems and shown to be able to achieve state-of-the-art results on complex combinatorial problems. Later, ant-based techniques were combined with the idea of using mobile agents for path management in next-generation telecommunication networks. The motivation was to apply stochastic optimization in the network and develop decentralized, adaptive, robust and self-organizing solutions to cope with the growing complexity of networks. On the basis of the promising results obtained for wired networks in the second half of the 1990's, ant-based systems for routing in multi-hop wireless networks (MHWNs) started to be developed, the first of them being ARA [START_REF] Günes | ARA-the ant-colony based routing algorithm for MANETs[END_REF] proposed in 2002.

In the context, ants are control packets that are used to repeatedly sample paths between source and destination pairs. At each node, the pheromone trails reflect the knowledge acquired by the colony, and the pheromone trail value associated with a given neighbour indicates the relative goodness of this neighbour to reach a specific destination.

The objective with this paper is neither to present yet another extensive survey nor a detailed description of the inner working of ant-based systems for MHWNs. For that, the interested reader is referred to [START_REF] Ducatelle | Principles and applications of swarm intelligence for adaptive routing in telecommunications networks[END_REF][START_REF] Farooq | Routing protocols for next-generation networks inspired by collective behaviors of insect societies: An overview[END_REF]. Rather, it is to look back on 10 years of research on applying ant-based techniques to routing in MHWNs and present some reflections on the challenges, evolutions, contributions, trends and perspectives in the field, touching upon both technical and non-technical aspects.

The rest of this paper is organized as follows. First, Section 2 addresses the application of ant-based techniques to mobile ad-hoc networks (MANETs), which is the most common type of MHWNs ant-based techniques have been applied to. Next, Section 3 discusses the evolution and current status of the research in the field and the perspectives for future developments and applications. Finally, concluding remarks are given in Section 4.

Applying Ant-Based Techniques to MANETs

This section addresses the application of ant-based techniques to MANETs. The focus is on unicast routing protocols as this represents the major part of the contributions, but there have also been proposals for using ant-based techniques for establishing and maintaining multicast trees in the context of MANETs, e.g. MANSI [START_REF] Shen | Ad hoc multicast routing algorithm with swarm intelligence[END_REF].

This section elaborates first on the premises for developing ant-based systems for routing in MANETs and the challenges related to the application of ant-based techniques in this context. The main solutions found in the literature are then presented. Finally, all these elements are summarized and discussed.

Premises

MANETs are autonomous MHWNs in which nodes communicate without the support of any pre-existing infrastructure and with no centralized control. Nodes may operate both as routers and hosts. They are all potentially mobile and battery-powered, and they may have limited processing and memory capacity. MANETs typically exhibit a high degree of dynamicity as nodes may move and leave or join the network at any time, making control, and in particular routing, a challenging task. Moreover, in addition to being able to deal with the changes in the environment, any control system designed for MANETs has to cope with the limited bandwith and the inherent unreliability of wireless communications, making this task even more challenging.

The main reasons that are put forward in the literature to support the application of ant-based techniques to routing in MANETs are: (i) the good match between the desirable properties of a path management systems for MANETs and the characteristics of ant-based systems in terms of adaptivity and robustness, and, as mentioned in the introduction, (ii) the promising results obtained applying ant-based techniques to routing in wired networks, in particular with AntNet [START_REF] Caro | AntNet: Distributed Stigmergetic Control for Communications Networks[END_REF]. Adaptivity is the ability of the system to adapt to changes and prevent disruptive events. This feature is empowered in ant-based systems by the repeated stochastic sampling of end-to-end paths and the resulting update of pheromone trail values at each node. It ensures the availability of timely routing information and enables the discovery and maintenance of multiple paths between source and destination pairs. Robustness is the ability of the system to maintain persitent operation under the occurrence of perturbations. Typical perturbations in MANETs include the loss of packets and the loss of routing information stored at nodes. Ant-based systems are robust in that they are not sensitive to the loss of individual ants and that routing information (pheromone trail values) is distributed throughout the network.

Challenges

The main challenge when applying ant-based techniques to routing in MANETs is efficiency. This is a two-sided challenge. Efficiency can be expressed in terms of routing overhead and in terms of convergence time. Others challenges include the proactive nature of ant-based systems, spreading the load, demands in terms of memory and processing and energy consumption.

Routing overhead

The routing overhead for ant-based systems is related to the repeated probing of paths realized by ants. This implies forwarding control packets back and forth between source and destination pairs which does consume bandwidth. As the wireless medium is shared and the available bandwidth is limited, controlling the overhead is therefore crucial for the performance of the system (data delivery ratio, delay, delay jitter). However, there is a trade-off between performance on the one side, and adaptivity and robustness on the other side. Ant-based systems are intrinsically robust to the loss of individual ants because they generate many of them. Reducing the number of ants directly affects the robustness of the system. In addition, the ability of the system to discover new paths and adapt to topological changes is conditioned by the fact that ants sample paths and update pheromone tables accordingly. Reducing the number of ants also directly affects the adaptivity of the system.

Convergence time Ant-based systems apply stochastic optimization and several iterations are required before such a system converges to a solution. Iterations corresponds to updates made by ants. Hence, the time it takes for a system to converge is related to the rate at which ants are generated. To be able to adapt, the time scale at which changes take place in the network should be much larger than that at which pheromone trail values are updated. Since MANETs are assumed to be highly dynamic environments and the number of ants required to keep up with the changes may be prohibitive with respect to the routing overhead incurred.

Proactive nature Ant-based systems rely on the proactive sampling of paths and comsume resources forwarding ants and maintaining pheromone tables at each node. In a wired networks, this makes sense as the topology is seldom changing, and, although nodes or links might failed, it is expected that they are restored with the same characteristics so that previously gathered data are still relevant and can actually be used to re-converge faster. In the context of MANETs, spending resources for proactive planning does generally not pay off as the conditions, in terms of topology, wireless link characteristics and traffic demands, are constantly changing. In addition, keeping historical data may be harmful for performance. When a node loses connectivity with one of its neighbor, although connectivity might be re-established at a later point in time with the same node, there is no guarantee that the quality of the path through this neighbor will be similar. Previously collected data may no longer be relevant and may lead to reduced performance if used.

Load spreading

The promising results obtained with AntNet in wired networks are in part due to its ability to achieve near-optimal load distribution by spreading the load along multiple path proportionally to their estimated quality. Achieving similar results in the context of MANETs is challenging for several reasons: (i) efficiently spreading the load in a highly dynamic environment is difficult, (ii) the behavior of the nodes and the traffic patterns in MANETs are unpredictable hence no assumptions that would help in distributing the load can be made, and (iii) distributing the load across multiple paths in MHWNs may lead to degraded performance because of interferences if the paths are not radio-disjoint [START_REF] Ducatelle | An analysis of the different components of the AntHocNet routing algorithm[END_REF][START_REF] Paquereau | Ant-based multipath routing for wireless mesh networks[END_REF].

Memory, processing and energy constraints Ant-based systems rely on the repeated and stochastic sampling of paths by ants and the maintenance of pheromone tables locally at each node. At each node, the pheromone table potentially contains an entry for each destination in the network through each of the neighbouring nodes. This implies repeated computations on a per-packet basis and imposes requirements in terms of processing and storage at each node. In addition, performing these operations and regularly forwarding ants are energy consuming tasks.

Solutions

Ant-based systems developed for MANETs are mostly derived from AntNet or directly apply the Ant Colony Optimization (ACO) meta-heuristic [START_REF] Dorigo | Ant algorithms for discrete optimization[END_REF]. Systems vary among others in terms of formulas used for pheromone updates, cost functions (hop-count, delay, or some measure of link quality including signalto-noise ratio, congestion or packet loss ratio, or any combination of those) and mechanisms used to improve the performance. Studies such as [START_REF] Baras | A probabilistic emergent routing algorithm for mobile ad hoc networks[END_REF][START_REF] Dhillon | Ant routing in mobile ad hoc networks[END_REF] have shown that directly applying ant-based systems developed for wired networks such as AntNet results in lower performance compared to conventional MANETs protocols, and in particular to AODV [START_REF] Perkins | Ad-hoc on-demand distance vector routing[END_REF], which is used as the reference protocols in most of the studies. In the following, the focus is on the mechanisms introduced specifically to address the challenges listed in Sect. 2.2.

Broadcasting forward ants Most of the systems broadcast forward ants to explore multiple routes concurrently and, hence, reduce the convergence time. However, this approach has also downsides: (i) paths are sampled regardless of their estimated quality, (ii) it results in a large overhead as ants get duplicated at each node, (iii) it results in the establishment of paths with low quality, and (iv) it leads to the discovery and reinforcement of a set of braided (interfering) paths that cannot therefore not be used to improve the performance of system by distributing the load. The resulting process is similar to the route discovery process of AODV, and some systems also introduce mechanisms used by AODV to limit the overhead such as expanding ring-search. To limit the overhead other systems such as AntHocNet [START_REF] Ducatelle | Adaptive Routing in Ad Hoc Wireless Multi-hop Networks[END_REF] resolve to broadcast only when no routing information is available at a node and otherwise apply a random proportional rule and unicast the ant only to one of the neighboring nodes.

Adopting a greedy deterministic forwarding policy Most of the systems abandon the stochastic forwarding policy used in AntNet for distributing the load accross multiple path and instead adopt a greedy deterministic forwarding policy, i.e. each node forwards data packets to its neighbor with the highest pheromone concentration for the given destination. The main reasons, as mentioned above, are that it is difficult to distribute the load in a dynamic environment and that distributing the load might lead to degraded performance because of interferences. Authors in [START_REF] Ducatelle | An analysis of the different components of the AntHocNet routing algorithm[END_REF] also show that the best performance with AntHocNet is achieved when a greedy deterministic policy is applied to the forwarding of ants. The reason is that the number of ants that can be generated is often too low for the system to be able to cope with the constantly changing environment, and the system cannot effectively explore multiple paths.

Relying on data packets for updating pheromone values In order to reduce the overhead, some of the proposed systems rely partly or solely on data packets to update pheromones, e.g. ARA and SARA [START_REF] Correia | Simple ant routing algorithm strategies for a (multipurpose) MANET model[END_REF]. The major drawback of this approach is that it affects directly the adaptivity and robustness of the system as it becomes dependent on data traffic to maintain timely routing information. Using data packets for exploration also results in decreased performance in terms for instance of increased delay and jitter.

Combining reactive and proactive path sampling

To reduce the overhead, most of the proposed systems adopt a reactive strategy, meaning that a node starts generating ants only when it is itself source of data traffic. Many systems keep, however, a proactive component and are therefore classified as hybrid, e.g. AntHocNet. The information collected by proactive ants may then used for bootstraping the reactive route discovery process and/or for improving the path and adapting to changes while data sessions are active.

Introducing additional control packet for link error notification and handling To improve the adaptivity of the system, many systems introduce additional dedicated control packets for explicit link failure notifications, local recovery and/or to deal with lost link failure notifications, e.g. ARA and An-tHocNet. One disadvantage of these mechanisms is that they create an extra overhead upon link failure. Another aspect is that local recovery techniques are opposite with the objective of converging towards a global optimum by constructing complete solutions, which is one of the advantages of ant-based systems as a stochastic optimization method.

Introducing pheromone diffusion This mechanism is introduced in AntHoc-Net. The idea is to let nodes share routing information with their 1-hop neighborhood via asynchronous beacon messages and derive virtual pheromones that are used to guide the exploration of ants.

Discussion

To deal with the challenges listed in Sect. 2.2 and bear the comparison with state-of-the-art reactive MANET routing protocols, in particular AODV, most of the proposed systems include mechanisms borrowed from those same protocols. This introduces a deviation from key elements of the optimization method (proactive sampling of end-to-end solutions, stochastic forwarding policy) and from what made its successful application to dynamic routing problem in wired networks (load-spreading). As a result, most of the systems are intermediate between conventional reactive MANET routing protocols and ant-based systems developed for path management in wired networks, such as AntNet. Some of them are much closer to AODV although the terminology used is that of antbased systems.

We argue that ant-based techniques have been applied to MANETs on wrong premises and question the suitability of this approach for routing in MANETs in the first place. The adaptability and robustness properties of an ant-based systems stem from the repeated sampling of path by ants, but are conditioned by the fact that the intensity of changes the system has to adapt to is sufficiently low compared to the rate at which ants sample paths and update pheromone trails. In addition, there is a trade-off between how adaptive and robust an antbased system can be and the routing overhead it generates. In typically settings for MANETs, changes occur too fast for the overhead not to be prohibitive.

Evolution, Status and Perspectives

In this section, we first present the evolution and the current status of the research on ant-based techniques for wireless networks. We then discuss when this approach is viable and when it is no, the question of the future of this approach and sketch possible future directions.

Evolution

The development of the first generation of ant-based systems for MHWNs, targeting MANETs, started at the beginning of the 2000's; the first ant-based system for MANETs reported in the literature being ARA [START_REF] Günes | ARA-the ant-colony based routing algorithm for MANETs[END_REF]. The explanation of why it started at that time is that: (i) routing in MANETs was a hot topic in the second half of the 1990's and early in 2000, and that, concurrently, (ii) promising results achieved by ant-based systems for wired networks were published by artificial intelligence researchers and attracted the attention of the MANET research community. Reversely, the development of ant-based systems for MANETs by networking researchers later attracted the attention of the artificial intelligence community (AntHocNet was first reported in 2004). Another reason is (iii) the observation made by De Couto et al. in 2002 [START_REF] De Couto | Performance of multihop wireless networks: Shortest path is not enough[END_REF] that minimizing the number of hops to decide on the best route between nodes in a MHWN is not necessarily synonym of maximizing the performance as longer links are generally more unstable than shorther ones. The fact that it is possible to integrate the link/path quality in the computation of the pheromone trail values is listed as an argument for applying ant-based techniques to MANETs in [START_REF] Günes | ARA-the ant-colony based routing algorithm for MANETs[END_REF].

Starting from 2004, and following the trend observed in MHWN research, a second generation of ant-based systems, targeting specialized MHWNs, started to be developed. The reason for this specialization is that the definition of MANETs is very generic and encompasses a wide range of scenarios and applications, which cannot all be best solved by a single protocol design. In particular, ant-based systems were developed for Wireless Sensor Networks (WSNs) [START_REF] Akyildiz | Wireless sensor networks: a survey[END_REF] and backbone Wireless Mesh Networks (WMNs) [START_REF] Akyildiz | Wireless mesh networks: a survey[END_REF]. The characteristics of these specialized MHWNs and of generic MANETs are summarized in Table 1.

- ferences MANET X X X X X WSN X X X X WMN X X X
In particular, a relatively large number of ant-based systems targeting WSNs have been designed, see for instance [START_REF] Saleem | A review of swarm intelligence based routing protocols for wireless sensor networks[END_REF] for a survey. WSNs have characteristics that make proactive approaches and load spreading a priori more appropriate than in generic MANETs: (i) nodes are static and, (ii) certain traffic patterns are expected between sensor nodes and sink nodes making load distribution both easier and more relevant. Compared to generic MANETs, restrictions in terms of memory and processing capacity as well as energy constraints are however more compelling. Hence, one can also question the intrinsic suitability of ant-based approaches in this context.

Similarly to WSNs, nodes in WMNs are static and the traffic is expected to be concentrated on certain paths (to and from gateway nodes providing connection to the Internet). On the other hand, contrary to WSNs, nodes in WMNs are gridpowered wireless routers with high processing and memory capacity. Therefore, although the number of contributions addressing specifically backbone WMNs is still limited, they seem to be the most actual type of MHWNs for the application of ant-based techniques. An example of such a system is for instance AMIRA [START_REF] Bokhari | AMIRA: interference-aware routing using ant colony optimization in wireless mesh networks[END_REF].

Finally, it should be noted that the research in the field has mainly been be driven by method-oriented rather than problem-oriented researchers. As for wired networks with AntNet, the most significant work in the field to date, AntHocNet, has been contributed by authors with a background in artifical intelligence. Applying ant-based techniques to MHWNs was, in a way, the next step after the promising results achieved with AntNet.

Status

As we assess the state of art, the research in ant-based techniques for routing in MANETs is about to cease, because: (i) as concluded in Sect. 2, ant-based systems are intrinsically not suitable for MANETs, and (ii) although several of the proposed systems and in particular AntHocNet are shown to outperform AODV, the gain in performance is not worth the added complexity. This latter assertion is supported by the work of AntHocNet authors who first observed in [START_REF] Ducatelle | An analysis of the different components of the AntHocNet routing algorithm[END_REF] that stochastically sampling paths results in lower performance compared to deterministically resampling the best paths, and who later proposed a new approach for integrating proactive and reactive routing in MANETs that built up on observations they made with AntHocNet, but with no references to ant-based techniques [START_REF] Ducatelle | A new approach for integrating proactive and reactive routing in MANETs[END_REF].

The application of ant-based techniques to specialized MHWNs is, on the other hand, still an open research topic. In spite of the fact that one may question its suitability for WSNs, many such systems have been proposed. It is, however, hard to assess the performance these proposals, since comparisons with state-of-the-art algorithms specifically designed for WSNs are seldom provided. Furthermore, both real-world implementations and testing, as well as mathematical modeling, is also listed as missing in [START_REF] Saleem | A review of swarm intelligence based routing protocols for wireless sensor networks[END_REF].

In the context of WMNs, few ant-based systems have been proposed so far. Pursuing these further, a first necessary step is a comparison with conventional state-of-the-art proactive routing protocols such as OLSR [START_REF] Jacquet | Optimized link state routing protocol for ad hoc networks[END_REF]. In the context of MANETs, ant-based systems have most of the time only been compared with AODV. To the best of our knowledge, no comparisons with traditional proactive MANET routing protocols have been reported.

Perspectives

Have ant-based techniques in MHWNs a future? If we expect it to be a simple, efficient and generally applicable technique, the answer is no. The technological rationale for this is found in Section 2. The method has hereto not been able to overcome the significant challenges put by the highly dynamical, wireless domain. On the other hand, we argue that the path management problem in WMNs is closer to that addressed in wired network than that addressed in MANETs, and that thus ant-based systems may have a future in this context. Moreover, there are specific aspects in the path management problem in WMNs for which antbased systems are relevant and attractive candidates. For instance, mesh routers are usually assumed to be equiped with multiple wireless interfaces and the channel assignment problem presents some similarities with the graph coloring problem for which good results have been reported using ant-based systems in a static context [START_REF] Costa | Ants can colour graphs[END_REF].

Although the path management problems in WMNs is closer to those of wired network than those addressed in MANETs, the challenges related to wireless communication -limited bandwidth and interferences -remain. Hence, applying systems developed for wired networks directly to WMNs leads to poor performance. One observation is that to rectify the problems incurred in the wireless domain, the typical approach has been to fall back on conventional path finding methods and patch ant-based systems with extra mechanisms, rather than devising new approaches. As an alternative direction for future research, we suggest to instead adapt the execution of the ant-based primitives to the context of wireless networks. One such approach is OAS [START_REF] Paquereau | Opportunistic ant-based path management for wireless mesh networks[END_REF]. The idea is to apply the opportunitic forwarding paradigm [START_REF] Bruno | Survey on diversity-based routing in wireless mesh networks: Challenges and solutions[END_REF] to the forwarding of ant agents. The fundamental changes are to defer the stochastic forwarding decision until after the reception of an ant by the potential next-hop nodes, and to execute it in a distributed manner among the receivers and in the time-domain rather than at the sender node and in the space-domain. This yields a reduction of the overhead and improves the adaptivity of the system compared to its unicast counterpart. Preliminary experiments indicate that such a system is able to achieve to efficiently spread the load along multiple paths [START_REF] Paquereau | Ant-based multipath routing for wireless mesh networks[END_REF].

Concluding Remarks

In this paper, we have retrospectively looked at 10 years of research in the field of ant-based systems for wireless networks. We have questionned the suitability of the approach for generic MANETs as well as for WSNs. We have presented the evolution and the current status of the research in this field and sketched a possible future direction of research for WMNs.

The suitability issue may also be raised for ant-based path management of dynamic telecommunication networks in general. In our view, the future might be in what we term ant-assisted systems rather than in ant-based systems. Rather than being the sole mean for self-organization, ant-based techniques seem more suitable as an element of a coumpound system using an array of techniques. In such a system, ants are not used for the primary path management task, but instead as a background process for the optimization of the primary system in the long run, e.g. the tuning of internal parameters, the detection of failure situations or degraded performance, or the optimization of external parameters affecting the path management problem. Applying ant-based techniques to the long-term optimization of the channel assignment in WMNs without relying on ants to solve the routing problem itself would for instance fall into this category.

Table 1 .

 1 Characteristics of MHWNs

	MHWNs	Limited	Mobility Predictable	High rate	Limited	Inter
		power,	traffic	of changes	bandwidth
		memory,	patterns	in the
		processing		network

. http://www.q2s.ntnu.no


Dharani Punithan
email: punithan.dharani@gmail.com

Bob Mckay

Self-organizing Spatio-temporal Pattern Formation in Two-dimensional Daisyworld

 was devised to demonstrate how the biota of a world could stabilise it, driving it to a temperature regime that favoured survival of the biota. The subsequent studies have focused on the behaviour of daisyworld in various fields. This study looks at the emergent patterns that arise in 2D daisyworlds at different parameter settings, demonstrating that a wide range of patterns can be observed. Selecting from an immense range of tested parameter settings, we present the emergence of complex patterns, Turinglike structures, cyclic patterns, random patterns and uniform dispersed patterns, corresponding to different kinds of possible worlds. The emergence of such complex behaviours from a simple, abstract model serve to illuminate the complex mosaic of patterns that we observe in real-world biosystems.

Introduction

Self-organization is a property of processes, in which patterns emerge at global scale due to the short-range local interactions of the system, with no external intervention [START_REF] Camazine | Selforganization in biological systems, 2nd Edition[END_REF]. A wide range of pattern formation processes have been observed in fields from chemistry (Turing's reaction-diffusion system [START_REF] Turing | The chemical basis of morphogenesis[END_REF], Belousov-Zhabotinsky reaction [START_REF] Belousov | A periodic reaction and its mechanism[END_REF]) and fluid dynamics (Bénard convection cells [START_REF] Bénard | Les tourbillons cellulaires dans une nappe liquide[END_REF]) to biology (patterns on animals' coats [6]).

In this paper, we investigate the emergence of self-organized, spatio-temporal vegetation patterns in two-dimensional (2D) Daisyworld [START_REF] Watson | Biological homeostasis of the global environment: The parable of daisyworld[END_REF]. First of all, we show that the surface temperature of planet is self-regulated around the mean 295.5K by the daisyworld dynamics of our model as in previous 2D daisyworld models [7,[START_REF] Ackland | Catastrophic desert formation in daisyworld[END_REF]. Then we observe the self-organizing properties of daisyworld, and qualitatively analyse pattern formation in daisyworld with Laplacian diffusion. We also quantify the emergence of vegetation patterns using a statistical measure, Moran's I. Pattern formation in daisyworld has been previously touched on in both 1D [START_REF] Adams | One-dimensional daisyworld: Spatial interactions and pattern formation[END_REF] and 2D [START_REF] Ackland | Emergent patterns in space and time from daisyworld: a simple coupled biosphere-climate model[END_REF] cellular automaton models. However our daisyworld model is extended to a coupled 2D chaotic system; we investigate daisyworld in diffusively coupled logistic systems in a 2D tordoidal regular lattice with chaotic local dynamics, a setting which has not been studied previously.

In our daisyworld model, depending on model parameter values we observe a variety of behaviours, such as 1) complex patterns (stationary/dynamic), 2) smoothly dynamic spatio-temporal patterns (Turing-like structures), 3) cyclic patterns (periodic spirals), 4) random patterns (stationary/dynamic) and 5) dispersed patterns (uniform).

Background

Daisyworld

Daisyworld was proposed by Watson and Lovelock [START_REF] Watson | Biological homeostasis of the global environment: The parable of daisyworld[END_REF] as a means to study planetary ecodynamics. This simple self adaptive system has one environmental variable (temperature) and two types of life (black and white daisies). The colour of the daisies influences the albedo (i.e. the reflectivity for sunlight) of the planet, and thereby influences the temperature of the planet. Though the daisies do not interact with each other directly, they do interact via the planetary temperature (environment); the behaviour of one type of daisies modifies the temperature, which affects the behaviour of the other, and vice versa. In general, increased growth of a particular type of daisy alters the temperature in a direction which is unfavourable for itself, but suitable to the other type. When black daisies start spreading more, the local temperature rises. Conversely, when white daisies spread, the temperature drops. Due to the contrary behaviours of the species, the temperature is self-regulated and daisies persist on the planet.

Turing instabilities

Turing instability of a homogeneous steady state in a two-species reaction-diffusion system provides theoretical mechanisms for pattern formation. Generally diffusion is perceived as an homogenising process. By contrast, Turing, in his remarkable paper (1952 [3]), demonstrated that a simple system of coupled reaction-diffusion equations could give rise to spatial patterns in chemical concentrations, through a process of chemical instability. The chemicals (morphogens or species), characterised as activator and inhibitor, react and diffuse throughout the tissue, and due to this physico-chemical process patterns such as spots and stripes can form in an animal's coat. These spatial patterns emerge from a homogeneous equilibrium state, due to symmetry-breaking via diffusion. This diffusion-driven instability forms the basis of pattern formation in physics, chemistry and biology. A wide spectrum of complex self organized patterns such as stationary, waves, spirals or turbulence emerge from reaction-diffusion systems.

Moran's I

Moran's I [START_REF] Moran | Notes on continuous stochastic phenomena[END_REF] is a measure of spatial autocorrelation, and is often used to analyse spatial patterns. The specific values -1, 0, +1 indicate spatially dispersed, random, and clustered patterns respectively. It is defined by equation 1:

I = C i j w ij i j w ij (P i -P )(P j -P ) i (P i -P ) 2 (1)
where C is the number of spatial units(cells) indexed by i and j, P is the size of the population, P is the mean population and w ij is a spatial weight matrix. In our work, w ij = 1 if i and j are horizontal or vertical neighbours, and w ij = 0 otherwise.

Self-organization and pattern formation in Daisyworld

Positive-negative feedback

Positive and negative feedback are the two basic modes of interactions among components of self-organizing systems [START_REF] Camazine | Selforganization in biological systems, 2nd Edition[END_REF]. Self-enhancing positive feedback coupled with antagonistic negative feedback leads to the formation of striking patterns in nature (ripples in sand dunes, schooling of fish, flocking of birds, coat of animals). A stationary state becomes unstable due to amplification, and this provokes pattern in the system. Since amplification purely in one direction leads to destruction, negative feedback takes control of the system and plays a critical role in inhibiting amplification and shaping the process and pattern, for example in colonial nesting of male bluegill sunfish [START_REF] Camazine | Selforganization in biological systems, 2nd Edition[END_REF].

In daisyworld, an initial growth of black daisies, due to positive feedback, increases the temperature of system and thus reinforces the change in the same direction. Hence it has the potential to explode the system. But in daisyworld, at a higher temperature, white daisies bloom and keep the positive feedback under control. In the same manner, the growth of white daisies, due to positive feedback, decreases the temperature in the same direction. But at a colder temperature, black daisies bloom and inhibit the system behaviour. The initial growth of daisies promotes new changes in the system state but on the other hand, increased growth counteracts the changes. The combined effect of positive and negative feedback couplings result in a negative feedback loop, which stabilises the physiological processes and leads to homeostasis [START_REF] Kump | The Earth System, 3rd Edition[END_REF]. Due to this behaviour, we can see that daisyworld is a self-organizing system.

Activator-inhibitor principle

Turing pattern formation in activator-inhibitor systems forms the paradigm of selforganization [START_REF] Nakao | Turing patterns in network-organized activator-inhibitor systems[END_REF], and provides a theoretical explanation of animal coat patterns [6]. The activator and inhibitor species destabilize the homogeneous state of the system due to diffusion, and lead to the spontaneous emergence of spatial patterns.

We can view daisyworld as a Turing reaction-diffusion system of two morphogens: black and white daisies. In cooler temperatures, black daisies absorb the available sun light and amplify their growth, functioning as activators (autocatalysis) because they activate their own production. But further increase in black daisies increases the surface temperature to the level at which cooler white daisies grow well, thereby acting as an inhibitor (self-suppression) because they inhibit their own production. This process is reversed for white daisies. Depending on the temperature, black and white daisies act both as inhibitor and as activator. Here, the diffusion of the two species, coupled by nonlinear local rules, result in the generation of patterns and spatial self-organization due to symmetry-breaking instability.

Model

The interacting components (cells) of the system are arranged on a toroidal 2D square (N × N ) regular lattice. Each cell is designated as a habitat with a maximum carrying capacity of 10,000 individuals. All sites in the lattice are randomly initialised with a population size in [0, 100] for both black and white daisies. The temperature is initialised randomly depending on the overlap between the optimal temperatures of black and white daisies (refer Table 2). The interaction rules among these system components are based only on local information (microscopic level) and not on the global emergence (macroscopic level). The black and white daisies, along with the temperature, diffuse around the world via their neighbourhood. The daisyworld model has two parts: the local component describing the effects on growth of daisies and change in temperature local to the particular cell; and the interaction component describing the migration of black and white daisies and diffusion of heat to their neighbours. Local Dynamics: The population size of each unit is computed based on the logistic equation with carrying capacity, refer equation 2:

P (t+1) = P t (1 + r[1 - P t K ]) (2)
where r is the bifurcation parameter (i.e. intrinsic capacity for increase). The parameter r provides positive feedback and the component [1-Pt K ] provides the negative feedback, with their combined effect regulating the population. Though the regulatory mechanism is built in, chaos emerges from the system [START_REF] May | Simple mathematical models with very complicated dynamics[END_REF]. The temperature change at each location is based on the energy balance equation [15] (refer equation 3):

σ B T 4 = SL(1 -A) (3)
where σ B is the Stefan-Boltzmann constant, T is the temperature, S is the solar constant, L is the luminosity and A is the albedo.

Global Emergence:

The behavioural rules in the system specify that black daisies grow better in colder temperatures, while white daisies grow better in hotter temperatures, both diffusing to their nearest neighbours. These simple components at the local level generate complicated dynamics at the global level due to the collective behaviour of the whole system, without the influence of any external forces. The interaction pattern is determined by the neighbourhood topology and the diffusion term. In our model, we use von Neumann neighbourhoods, consisting of a central cell and its four orthogonal neighbours (top, down, left and right, without diagonal interactions), and a Laplacian diffusion model for both species and temperature. The key ingredients leading to global emergence are the local nonlinear interactions among the components, physiological and behavioural rules, and physical constraints. The inhibition of self-growth can also arise from the physical constraints, due to the self-limiting factor of carrying capacity.

Albedo:

The albedo of the planet can be computed by equation 4:

A = A b α b + A w α w + A g α g (4)
where α b , α w , α g (= 1α wα b ) ∈ [0, 1] are the relative areas occupied by black, white daisies and bare ground, A b is the albedo of ground covered by black daisies, A w that of ground covered by white daisies and A g that of bare ground. We assume that A w > A g > A b , with corresponding values of 0.75, 0.5, 0.25.

Temperature:

The local temperature update is defined by diffusion, heat radiation, solar absorption and Gaussian noise, as in equation 5:

C δT (i,j,t) = D T ▽ 2 T (i,j,t) -σ B T 4 (i,j,t) + SL(1 -A (i,j,t) ) + Cǫ (i,j,t) (5)
where C = 2500 is the heat capacity, δT i,j,t = T i,j,t+1 -T i,j,t is the change in temperature at (discrete) time t, 1 ≤ (i, j) ≤ N indicates a 2D lattice point, N is the lattice length, D T = 500 is the diffusion constant, ▽ 2 T is the Laplacian operator, σ B is the Stefan-Boltzmann constant, S is the solar constant, L is the luminosity, A is the albedo and ǫ represents Gaussian white noise (with mean zero and standard deviation 1.0) multiplied by the noise level.

Growth:

The growth rate of daisies is defined as a parabolic function as in equation 6:

β(T ) = 1 -[ (T opt -T ) 2 17.5 2 ] ( 6
)
where T is the local temperature, and T opt is the optimal temperature of the species, which depends on their colour. The optimal temperature for black daisies is lower than for white. In general, growth rates of daisies depend on local temperature, which in turn depends on the albedo and thus on the proportions of the daisy species and bare ground.

Population Size:

The local population update is computed based on Laplacian diffusion and local population growth, governed by the logistic equation with carrying capacity, as in equation 7:

δP (i,j,t) = D ▽ 2 P (i,j,t) + rP (i,j,t) [β(T ) - P (i,j,t) K ] (7)
where at location (i, j) and time t, P i,j,t is the population size, D is the fraction of the population being dispersed to its neighbours, r is the natural rate of increase, β(T ) is feedback coefficient, and K is the carrying capacity. Deterministic chaos can also be induced by diffusion and hence our system doesn't require high intrinsic growth rates.

Assumptions and Limitations

The model parameters are represented in Table 1. They reflect the general daisyworld literature and parameter values [START_REF] Watson | Biological homeostasis of the global environment: The parable of daisyworld[END_REF]7]. The last two rows of Table 1 are variable parameters, while the rest are fixed parameters. The range of temperatures in which daisies can survive is the viability constraint of our model, and is based on the overlap between the optimal temperatures of black and white daisies as in Table 2. We define the optimal temperature of black daisies within the range [278, 295.5]K and of white in the range [313, 295.5]K so as to have an overlap in the range [0, 100]%. For this work, we use a 10% overlap for most of experiments, with a 0% overlap for one. As we have a huge parameter space, we have restricted our experiments to an increase rate of r = 1. We used the same diffusion constant for both black and white daisies in all our experiments. We limited our focus to studying the effects of species diffusion (which influences the species distribution) and noise (which influences the temperature) in detail.

Results

The results presented here are samples from over 15,000 experiments with different variable parameter settings, including species increase rates in the range [START_REF] Watson | Biological homeostasis of the global environment: The parable of daisyworld[END_REF][START_REF] Belousov | A periodic reaction and its mechanism[END_REF], species dispersion rates in {0, 0.001, 0.01, 0.1, 0.2, 0.5, 1}, overlap in {0, 5, 10, 25, 50, 100}%, noise level in {0, 0.001, 0.005, 0.01, 0.05, 0.1, 0.2, 0.5, 1, 2.5, 5, 7.5, 10}, and temperature diffusion in {0, 0.05, 0.1, 0.2, 0.3, 0.5, 1}. We examined the results for the emergence of interesting patterns both qualitatively (by capturing snapshots of each experiment and inspecting visually) and quantitatively (using the value of Moran's I). Our preliminary investigation found interesting behaviours in the 0% and 10% overlap settings, and we present these here. Our future aim is to develop suitable metrics to automatically select interesting cases from the parameter space. We present the results by showing the state of the Daisyworld at different epochs. Since the behaviour is different in each case, the epochs chosen differ in each case. In these visualisations, a location is shown as black if the population of black daisies at a particular location is larger than that of white, and vice versa. The variable parameter values used in each case are detailed in the caption of the figure.

As it is impossible to show all spatial snapshots over 5000 epochs, we have plotted the temporal dynamics of the global temperature, which characterises the nature of the daisyworld. We note that in all cases, the temperature self-regulates around 295.5K, confirming the persistence of the daisies as an emergent property of daisyworld. We have also plotted the value of Moran's I for both black and white daisies.

A wide spectrum of patterns emerge from our model, depending on the model parameter values. We have classified them into 1) complex patterns (stationary/dynamic), 2) smoothly dynamic spatio-temporal patterns (Turing-like structures), 3) cyclic patterns (periodic spirals), 4) random patterns (stationary/dynamic) and 5) dispersed patterns (uniform). We were able to find examples of most behaviours with an overlap of 10%; however for Turing-like structures, we used an overlap of 0%. A maze-like self-organizing pattern (refer Fig. 1) emerge from our model with a diffusion constant 0.001 for both species and a noise level of 0.001. Though the ini-tial states look chaotic, the system finally converged to a stationary maze-like structure. The global spatial behaviour is plotted in Figure 2. The temperature initially fluctuates, but then converges, exhibiting a stationary pattern formation. Initially Moran's I is very high, indicating the formation of very large clusters, but it finally drops below 0, indicating a uniformly dispersed pattern. In between, the fluctuations show complex dynamics. Raising the diffusion constant to a still relatively low 0.1 for both species, and with a noise level of 0.001, we see the emergence of a 'modern art' complex pattern (Figure 3); Figure 4 shows its global dynamics. Moran's I stays very high throughout the run till the final epoch -and visually we can see clusters in some locations in that final epoch.

Self-organized complex pattern

Turing-like structures

Turing-like patterns emerge with a diffusion constant of 0.5 for both species, a noise level of 0.2 and an overlap of 0% -refer to Periodic spirals can be seen in the behaviour when we use a diffusion constant of 0.2 for both species and a noise level of 0.05 (Figure 7). The corresponding global dynamics are plotted in Figure 8. The cycles in temperature seen in sub-figure (a) confirm the periodicity in pattern formation. Except for the initial states, Moran's I is high throughout the run, confirming the formation of large clusters (sub-figures (b) and (c)).

Random pattern

Random patterns evolve from the system for a noise level of 0.01 and with no species dispersion -see Figure 9 Interestingly, we do see the emergence of a hint of a large scale grid in this setting, presumably resulting from the effects of temperature diffusion.

Conclusions

Summary

The results illustrate that the underlying dynamics of daisyworld are not restricted to stationary patterns but extend to periodic and chaotic behaviour. We see stability as The localised interaction of species (black and white in our daisyworld model generate an array of fascinating spatio-temporal patterns, due to endogenous selforganization resulting from "symmetry-breaking instability". In our nonlinear reactiondiffusion system, symmetry-breaking is induced by species dispersion, temperature diffusion and noise. A low species dispersion rate combined with noise leads to complex self-organized patterns (Figures 1 and3). As species dispersion and noise increase, periodic (Figure 7) and dynamic chaotically evolving patterns (Figure 5) emerge. In the absence of species dispersion, random (Figure 9) and dispersed (Figure 11) patterns emerge. We thus see the constructive roles of spatial diffusion, nonlinear local rules and noisy fluctuations in influencing the physiological, behavioural rules and physical constraints and thereby self-organizing the daisyworld.

Future Work

Since our parameter space is huge, it is practically infeasible to categorise all domains of pattern formation in daisyworld by qualitative visualisation -some automated criterion is required. We used Moran's I, but it assists only with spatial analysis, omitting the temporal analysis. In order to to fully analyse both parameter and pattern formation space, we need objective statistical measures which can uniquely classify the patterns in both space and time. We hope that further work will help us to understand the intrinsic connection between spatio-temporal dynamics and pattern formation in daisyworld.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Maze like organization for D = 0.001 and noise level= 0.001 in 2D 100 × 100 Lattice

Fig. 3 .Fig. 4 .

 34 Fig. 3. Modern art complex pattern for D = 0.1 and noise level= 0.001

Figure 5 .

 5 The corresponding global dynamics are plotted in Figure6. The patterns which emerge are smooth and slowly changing, looking like an animal's coat (e.g. Holstein cow). The fluctuations in temperature in sub-figure (a) of Fig.6show the dynamicity of pattern formation. The high Moran's I throughout the run (except the initial transitions) corresponds to the emergence of large clusters in the patterns (see subfigures (b) and (c) of Figure6).

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Turing-like structures for D = 0.5, noise level= 0.2 and overlap0%

Fig. 8 .Fig. 9 .

 89 Fig. 8. Global dynamics of periodic spirals

Fig. 10 .Fig. 11 . 5 TimeFig. 12 .

 1011512 Fig. 10. Global dynamics of random pattern

Table 1 .

 1 Daisyworld Parameter Settings (last two rows are variable parameters, others are fixed)

	Parameter	Value(s) Parameter	Value(s)
	Fixed
	Number of cells	100 × 100 Bare ground Albedo (Ag)	0.5
	Heat Capacity (C)	2500 Albedo of black daisies (A b )	0.25
	Diffusion constant (DT )	500 Albedo of white daisies (Aw)	0.75
	Stefan-Boltzmann constant
	(σB) E -8 W m -2 K -4	5.67 Mean growth temperature	295.5 K
	Luminosity (L)	1 Solar Insolation (S) W m -2	864.65
	Variable
	Diffusion rate of black daisies	[0, 1] Diffusion rate of white daisies	[0, 1]
	Noise Level	[0, 10] Natural Rate of Increase	[1, 4]

Table 2 .

 2 Optimal temperature overlap

	Overlap(%) Opt. temp. (K) Opt. temp. (K) Initial values (K)
		black daisies	white daisies
	0	278	313	[270, 320]
	10	284.5	306.5	[280, 310]

Acknowledgements

This research was supported by the Basic Science Research Program of the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Project No. 2011-0004338), and the BK21-IT program of MEST. The ICT at Seoul National University provided research facilities for the study. We would like to thank Ilun Science and Technology Foundation for their generous support to Dharani Punithan.


Efrat Blumenfeld-Lieberthal
email: efratbl@post.tau.ac.il

Juval Portugali

MetroNet: A Metropolitan Simulation Model Based on Commuting Processes

Keywords: urban complexity, urban networks, agent based models, urban simulation models

The aim of this work is to identify a set of fundamental rules that govern the interactions within urban systems at the metropolitan scale. For that, we developed an USM (Urban Simulation Model) specifically designed to study the evolution and dynamics of systems of cities. Our model is innovative in its structure: it is a superposition of cellular automata and agent based modeling approaches (that are essentially spatial analyses) and a complex network approach (that is essentially a topological analysis). This implies that in our model, the local activities and interaction of agents give rise to the global urban structure and network that in turn affects the agents' cognition, behavior, movement and action in the city and so on in circular causality. The model simulates commuting patterns of agents within a metropolis. The agents in our model represent workers who look for working places, the nodes represent urban employment centers, and the links represent commuters. Our results address three issues: the first suggests that the perception of urban boundaries plays a significant role in the metropolitan evolution in terms of network topology. This means that the existence of business centers, located in proximity to each other (but belonging to different municipalities) may lead to the emergence of new centers at the metropolis scale. The second issue concerns urban segregation; our results suggest that the location preferences of the agents regarding proximity to similar/different agents have a major affect not only on the urban morphology but also on the topology of the urban network. The third and last issue concerns the size distributions of agents in our model; these distributions correspond to all types of homogenous distributions observed in real system of cities.

Introduction

In 1965, Christopher Alexander in his paper "A city is not a tree" described the differences between the way modernist planners grasp the relationships between urban entities and the way they exist in reality. A careful reading of this work suggests that Alexander was one of the first scholars that proposed a new understanding of cities as complex systems that can be described as complex networks. These networks contain many entities that interact in many levels and scales. According to Alexander, they must be treated as a system rather than separate objects that act within the urban context.

When looking into cities, it always comes down to urban agents that interact with each other. These agents can be individual human beings, families, households, organizations, municipalities, firms, and so on. They interact in a variety of domains such as economic, infrastructures, transportation facilities, commuting, internal immigration and more.

Since the works of Watts and Strogatz [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF] and of Barabasi and Albert [START_REF] Barabási | Emergence of Scaling in Random Networks[END_REF], the science of complex networks has been developing rapidly. It has been applied to various disciplines in order to study the topology of large networks and to understand their development and robustness [START_REF] Albert | Error and attack tolerance of complex networks[END_REF][START_REF] Gallos | Stability and topology of scale-free networks under attack and defense strategies[END_REF][START_REF] Paul | Optimization of network robustness to random breakdowns[END_REF][START_REF] Jeong | The largescale organization of metabolic networks[END_REF][START_REF] Newman | Random graph models of social networks[END_REF][START_REF] Jeong | Measuring preferential attachment in evolving networks[END_REF][START_REF] Blumenfeld-Lieberthal | The topology of transportation networks: a comparison between different economies[END_REF][START_REF] Patuelli | The evolution of the commuting network in Germany: spatial and connectivity patterns[END_REF][START_REF] Batty | Cities and complexity[END_REF]. Most of these works, however, consider the topological characteristics of urban networks, while the physical aspect is mostly neglected.

Cities are typical examples of complex, self-organizing systems [START_REF] Patuelli | The evolution of the commuting network in Germany: spatial and connectivity patterns[END_REF][START_REF] Batty | Cities and complexity[END_REF][START_REF] Porutgali | Complexity, Cognition and the City[END_REF]. They have originally emerged and are still developing out of the interactions between many agents that are located and move in space and time [START_REF] Portugali | Self-Organization and the City[END_REF][START_REF] Alexander | A city is not a tree[END_REF]. These agents are motivated by a variety of forces ranging from cognitive capabilities and needs to economic considerations, political ambitions, etc., with no central force that affects their behavior. These interactions entail a huge number of links that create complex networks which form the city.

Most if not all complexity theories and models have been applied to the study of cities with the implication that we now have a whole family of urban simulation models that model cities (See review in [START_REF] Batty | Cities and complexity[END_REF][START_REF] Batty | Fractal Cities: A Geometry of Form and Function[END_REF]). In the last decade or so cellular automata and agent-base models (CA/AB) have become the main medium to simulate cities as complex systems [START_REF] Batty | Fractal Cities: A Geometry of Form and Function[END_REF][START_REF] Benenson | Geosimulation: Automata-based modeling of urban phenomena[END_REF][START_REF] White | Cellular automata and fractal urban form: a cellular modeling approach to the evolution of urban land-use patterns[END_REF] while in the last few years we see studies that model cities as complex networks that are often typified by power law distributions [START_REF] Porutgali | Complexity, Cognition and the City[END_REF][START_REF] Andersson | A complex network approach to urban growth[END_REF][START_REF] Porta | The network analysis of urban streets: A dual approach[END_REF].

In this work, we focus on urban networks and study their topological as well as spatial characteristics. We aim to identify a set of fundamental rules that govern the interactions within urban systems at the metropolitan scale. For that, we developed an USM (Urban Simulation Model) specifically designed to study the evolution and dynamics of systems of cities. Our model is innovative in its structure as it is a superposition of cellular automata and agent based modeling approaches (that are essentially spatial analyses) and a complex network approach (that is essentially a topological analysis). This implies that in our model, the local activities and interaction of agents give rise to the global urban structure and network that in turn affects the agents' cognition, behavior, movement and action in the city and so on in circular causality.

In the next section we present the detailed description of the above model. Then, we present some preliminary results. As the presented work is an ongoing one, in the last section, we elaborate on directions for future work.

A detailed description of the model

The essence of our model can be described by the following scenario: we start with a metropolitan environment divided into sub-areas, some of which represent cities with their municipal boundaries and others represent green areas. Each city is growing logistically (e.g. demographically due to natural increase and/or migration), and economically due to the development of employment and business centers. The residents (agents) of the city seek working places either in their town of residence or in other cities (to which they need to commute). Their decision where to work in based on a gravitation spatial interaction model i.e. their choice of a preferred working place is proportional to the size of the existing employment centers and inversely proportional to the distance between the agent's home and employment center. The choices where to work give rise to an intra-and inter-city commuting system. When a certain threshold of commuters between two cities is crossed, a link between these cities is created. The weight of the link is a dynamic parameter, which represents the volume of commuters.

Our USM is composed of the following elements and process:

1. The metropolis infrastructure: (a) The area of the metropolis which is represented by a rectangular is defined (by the user) (b) The model divides the metropolis into 25 spatial equal units and characterizes 7 of them (based on a random choice) as green areas (where urban development is restricted) (c) The remaining 18 urban units are randomly united into 12 municipalities. These numbers are comparable to the Tel Aviv metropolis (to which we intend to relate this work in later stages). (d) For each municipality, a center is randomly selected and being occupied by one of the agents. This random choice of the urban center can be explained as a historical accident [START_REF] Porutgali | Complexity, Cognition and the City[END_REF]. 2. Characterizing up to 3 types of agents. Each type of agents is represented by: (e) Growth rate -the addition of new agents in each iteration (f) Sensitivity to size of employment opportunities (α -see below) (g) Sensitivity to distance between the agent's place of residence and work (β -see below) 3. Location method:

(h) At the first stage: for each municipality an agent is randomly selected and located. Each new agent can be located either in its original city or elsewhere based on the availability of vacant space and a probability function. If the agent is not located in its original town its new location is determined by a gravita-tion function. The gravitation function can be described by: ܩ = ሺெሻ ഀ ഁ , where G represents the gravitation force, M and m represent the size (total number of agents) in the agent's original and destination towns (correspondingly), r represents the physical distance between the centers of both towns, and α and β are parameters that characterize the sensitivity of the agent to size and distance correspondingly. (i) The segregation function can be turned On/Off. This determines whether agents want to locate close to agent of their type. 4. Network definition:

(j) Nodes are defined at the location of the first agent in each city. The historical center is considered as the city's core (k) A minimum volume of commuters (agents) is needed to define a link. This number is predefined as a percentage of the population (and is changeable) (l) The weight of the links is presented (gradient color change) (m) When a predefined percentage of the cells, located in proximity to the municipal boundaries are occupied by agents, an agglomeration emerges. This means that the municipal boundaries are ignored in the gravitation function and the new center of the agglomeration is moved to its geographical center. Note that this condition is ignored when the predefined percentage of cells, located near the municipal boundary is 100.

Preliminary results

Based on the above USM, we can introduce some interesting preliminary results. In figures 1-2 (bottom), the green rectangles represent green areas where urban development is restricted, while the black lines represent municipal boundaries. The different types of agents are represented by black, red, and grey pixels. The first issue we address concerns agglomerations in face of municipal boundaries. Figure 1 presents the morphology and the network topology of the metropolis in two different runs of the model. In both runs there is only one type of agents. In the first run of the model (figure 1a), the agents consider the municipal boundaries in their choice of working places, while in the second run (figure 1b) the agents consider agglomerations instead of municipalities in their choice of working places. It can be seen that the perception of urban boundaries plays a significant role in the metropolitan evolution in terms of network topology. While in figure 1a The second issue concerns urban residential segregation. Figure 2 presents the morphology and network topology of the metropolitan in two different runs of the model in which there are 3 different types of agents. In the first run the agents are indifferent to the kind of their neighbors, thus, they might be located in adjacency to other kinds of agents (figure 2a, bottom). In the second run the agents prefer to locate in proximity to agents of their own kind (figure 2b bottom). To control this behavior, we added a segregation function to the model. When in use, this function affects the gravitation function such that M and m are calculated based only on the number of agents of the same kind as the agent who looks for a location.

It can be seen that the location preferences of the agents regarding proximity to similar/different agents have a major affect both on the urban morphology and network topology. When the agents prefer to locate in proximity to their own kindurban segregation emerges. In terms of network topology, a significant difference between both runs is observed. When plotting the rank-size distribution of the links' weight (Y=log(weight of link) and X=log(rank of link)) the distribution of the first run (agents have no preferences of neighbors) corresponds to an exponential function (figure 3a), while the distribution of the second run (agents prefer to work in proximity to their own kind) obeys a power law (figure 3b).

In addition, in the first run the highest level of commuters reaches approximately 11% of the total population in the metropolis. On the other hand, in the second run the highest level of commuters reaches only 6.7% of the total population. This can be explained by the fact that in the first run, the large employment centers attracted all three types of agents. However, in the second run, employment centers attracted only agents of the same type as their majority population. Lastly, we address the rank-size distribution of the agents in our model. There is an extensive work on scaling relationship in cities [START_REF] Jiang | Topological analysis of urban street networks[END_REF][START_REF] Zipf | National unity and disunity[END_REF][START_REF] Krugman | The Self-Organizing Economy[END_REF][START_REF] Batty | Rank clocks[END_REF][START_REF] Geoffery | Synchronous neural activity in scale-free network models versus random network models[END_REF][START_REF] Benguigui | A new classification of city size distributions[END_REF]. One of the classic works is that of Haggett [START_REF] Haggett | Models in geography: the Madingley lectures for 1965[END_REF] who suggested that city size distributions of urban systems can be divided into three classes: the first includes system of cities with primate city (the largest city is considerably larger than the next largest one), the second class includes power low distributions and the third includes systems of cities in which the large cities are rather homogenously distributed. An empirical study on the size distribution of cities [START_REF] Benguigui | A new classification of city size distributions[END_REF] showed that homogenous systems of cities can indeed be classified into these classes. When analyzing the resulted size distributions of the different agents in our model, we get all the observed population size distributions (see examples in figure 4). This, of course corresponds to the different values of parameters we used in different runs. This is also valid to the network topology. We found that the rank size distributions of the links' weight also correspond to these three classes. These results suggest that the set of laws to which the agents in our model obey, correspond to the behavior of people in real systems of cities. By following these laws and the values of the parameters that represent them we can try and understand which forces are the most dominant in the creation of scaling relationship within urban systems (as opposed to other parameters that lead to other relations).

Discussion and conclusion

In this paper we've presented an urban simulation model designed to study the evolution and dynamics of metropolitan systems of cities. Our model is built as a superposition of cellular automata and agent based modeling approaches (that are essentially spatial analyses) and a complex network approach (that is essentially a topological analysis). We study both the topology and the spatial characteristics of urban commuting networks at the metropolis scale. In a future work, we intend to elaborate our model in order to study additional urban phenomena; the first addresses the morphology of the municipal boundaries. To explore this issue we will change the definition of municipal boundaries and study the effect of increased irregularity of the boundaries' morphology on the resulted networks. The second phenomenon addresses preservation of rural environment under the economic pressures of real-estates developers. In reality rural areas are often transformed into developed urban environments (e.g. edge cities). In our next version of the model, we will enable (under controlled conditions) urban development in the rural (green) areas and study its effect on the urban networks in terms of the networks' topology and spatial characteristics.

Fig. 1 .

 1 Fig. 1. The model at the metropolitan scale with one type of agents: a) municipal boundaries are considered and b) agglomerations of built-up areas replace municipalities

Fig. 2 .Fig. 3 .

 23 Fig. 2. The model at the metropolitan scale with 3 types of agents: a) agents are tolerant to location in proximity to other types of agents b) agents prefer to locate in proximity to their own types

Fig. 4 .

 4 Fig. 4. The three classes of city size distributions, recovered by different runs of our model (with different parameters at each run)