Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://chaos.phys.msu.ru/ryabov/pdfs/CoevMotion2008.pdf Äàòà èçìåíåíèÿ: Wed Jul 15 14:50:49 2009 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 19:56:43 2012 Êîäèðîâêà:

Eur. Phys. J. Special Topics 157, 223­238 (2008) c EDP Sciences, Springer-Verlag 2008 DOI: 10.1140/ep jst/e2008-00643-9

THE EUROPEAN PHYSICAL JOURNAL
SPECIAL TOPICS

Co evolutionary motion and swarming in a niche space mo del of ecological sp ecies interactions
C.J. Dommar1 , a , A. Ryabov1,2 , and B. Blasius2
1 2

,b

Department of Physics, University of Potsdam, 14415 Potsdam, Germany ICBM, University of Oldenburg, 26111 Oldenburg, Germany

Abstract. Organisms are involved in coevolutionary relationships with their competitors, predators, preys and parasites. In this context, we present a simple model for the co-evolution of species in a common niche space, where the fitness of each species is defined via the network of interactions with all other species. In our model, the sign and type of the pairwise interactions (being either beneficial, harmful or neutral) is given by a pre-determined community matrix, while the interaction strength depends on the niche-overlap, i.e. the pairwise distances between species in niche space. The evolutionary process drives the species toward the places with the higher local fitness along the fitness gradient. This gives rise to a dynamic fitness landscape, since the evolutionary motion of a single species can change the landscape of the others (known as the Red Queen Principle). In the simplest case of only two-species we observe either a convergence/divergence equilibrium or a coevolutionary arms race. For a larger number of species our analysis concentrates on an antisymmetric interaction matrix, where we observe a large range of dynamic behaviour, from oscillations, quasiperiodic to chaotic dynamics. In dependence of the value of a first integral of motion we observe either quasiperiodic motion around a central region in niche space or unbounded movement, characterised by chaotic scattering of species pairs. Finally, in a linear food-chain we observe complex swarming behaviour in which the swarm moves as a whole only if the chain consists of an even number of species. Our results could be an important contribution to evolutionary niche theory.

1 Intro duction
The concept of an ecological niche is one of the most discussed, rethought and reviewed sub jects in biology [1]. While most ecologists today consider niche theory as a central and important concept which serves as a link unifying evolutionary, population, community and ecosystem ecology, the controversial discussion about the definition of a niche still remains a hotspot in the scientific community ­ even to the extreme of some authors who propose that the concept may not be needed at all in order to explain ecological and biogeographical patterns [2]. There are two traditional schools defining the niche concept from different points of view: one referring to the habitat a species needs for survival, and the other referring to the role or impact the species plays in the community [1]. For the first time, the term "niche" was explicitly used by Grinell [3] who defined a niche as "all the sites where organisms of a species can live", i.e., where conditions are suitable for life. Ten year later, Elton [4] considered a niche as a function performed by a species in its community, specifically the species's functional role within the
a b

e-mail: carlos.dommar@gmail.com e-mail: blasius@icbm.de


224

The European Physical Journal Special Topics

Table 1. List of all possible types of two-species interaction as represented by the elements aij of the interaction matrix A. aij 0 -1 +1 -1 +1 +1 aji 0 0 0 -1 +1 -1 Type of interaction Neutralism, complete lack of interaction Amensalism, the interaction is detrimental to one species but neutral to the other Commensalism, one species benefits and the other is neither benefited nor harmed Competition, the interaction is harmful for both species Mutualism, both species experience benefit from each other Predation or parasitism, one species benefits and the other is harmed

food chain. Thus, the first definition emphasises the "address", while the second one stresses the "profession" of a species [5]. Later on, Hutchinson [6] extended the concept, putting it in the framework of mathematical formalism. His definition for a niche was that of the "sum of all environmental factors" acting on the organism, hence the niche space would be an "m -dimensional hypervolume". In this way Hutchinson defined any number m of limiting factors (e.g., temperature, height, resource, salinity, etc.) for a given organism and the niche is regarded as the range of those n factors in which the species can live. On the other hand, the Red Queen evolution hypothesis, as proposed by Van Valen [7], states that for a species in an ecosystem it is necessary to develop continuously just in order to maintain its fitness relative to the other evolving species. Thus the system can evolve without any external influence, but only due to internal biological factors. The more interaction among the species the stronger the Red Queen effect. This intriguing observation opened a polemic discussion in the evolutionary biology regarding to the modes of evolution, i.e. continuous evolution, gradualism and punctuated equilibrium [8, 9], the evolution of sex [10, 11] and ecological communities assembly [12]. The Red Queen notion has been investigated in the context of the evolutionary ecology [13, 14] and adaptive dynamics [21­23]. In this article, based on Hutchinson's metaphor of a hyper-dimensional niche space we investigate Red Queen evolutionary dynamics of a community of n species. We define the species as points in an m-dimensional niche space. The ecological relationships are defined by a pre-determined community matrix, while the interaction strength results from the pairwise distances of the species in the niche space. In this setting the coevolutionary dynamics is implemented by the maximisation of the fitness function. This gives rise to an intricate coevolutionary motion in niche space, which is investigated with analytical and numerical methods. We first investigate the simplest case of only two species. Then we explore all anti-symmetric three species communities and finally we study n-species predator-prey chains.

2 The mo del
Consider an m-dimensional niche space . The dimensions of may represent any environmental parameter or living condition of a species, e.g. temperature, humidity, amount of light, depth, etc. Our model describes the evolutionary dynamics of n species in this niche space, where each species i = 1 ...n is characterised solely by its position xi = {xi1 ,...,xim }. The ecological role of a species within the community is determined by the matrix A with elements aij {-1, 0, +1}, which define the influence of the i-th species on the j -th one. The matrix A thus describes the adjacency matrix of the complex network of species interactions and for example, might be thought to represent the topology of a food-web. The pairwise coefficients aij and aji complete all classical types of ecological interactions between two species (see Table 1). Note, that we do not consider self-interactions, i.e., the diagonal of A is set to zero, aii = 0. Further, in our model the community matrix A is fixed a priory and does not change in the course of time. While the sign or the presence of the interaction between any pair of species is characterised by the matrix aij , we assume that the interaction strength (xi , xj ) is a function of the species


Active Motion and Swarming

225

positions in . The main idea here is that two species are able to interact more strongly if they are close together in niche space, while they cannot interact if they occupy different niches. Thus the interaction strength should decay with the pairwise distance dij = |xi - xj | in niche-space, i.e. = (|xi - xj |), independent of the interaction type. Of course, not all species which are close in niche space do actually interact at all. Also, the proximity of two species in does not tell us anything about the type of interaction (e.g., predator-prey, competition, etc.). Such properties in our model are determined by the matrix A. Thus, we propose to separate the ecological interactions into two parts: a time-independent community matrix, which gives the signs and type of the interaction, and a time dependent function (|xi - xj |), which depends on the positions in niche space and determines the interaction strength. With this in mind we suggest that the evolutionary fitness Fi of a species i is composed as a sum of the pairwise contributions from all other species. Each contribution is given as a product aij (dij ) of the elements of the interaction matrix aij (the species' ecological role) with the distance dependent interaction strength (dij ) (the species's similarity)
n

Fi =
j =1

aij (|xj - xi |)+ fi (xi ),

i = 1,...,n.

(1)

Here fi (xi ) describes an additional, external contribution to the fitness of species i which specifies the favourable region or the niche preference in the absence of all other species. It could represent for example a radial cost. Note, that the sum in Eq. (1) runs over all species j = i. In other words, we preclude evolution of species self-interaction and set the diagonal elements to zero, aii = 0. A common suggestion is that the interaction strength between two species should be a monotonically decreasing function of the distance dij separating these species in . In particular, in this presentation we model it as a Gaussian with a common width for all species pairs -d2 ij (dij ) = exp . (2) 2 2 Renormalizing the niche space, we can set = 1 without loss of generality. Note that while the interaction strength is defined to be strictly positive, the fitness Fi of a species in Eq. (1) may well be negative. Having with Eq. (1) defined the fitness Fi of a species as a function of its position in , it is straightforward to specify the resulting evolutionary dynamics. We assume that each species evolves along the gradient of its fitness landscape trying to increase its adaptation within the community. Thus we propose the following dynamic equations xi = i i Fi , (3)

where i denotes the mutation rate of species i. Using Eq. (1) the dynamics of each species then is given by the set of evolution equations 1 dxi = i dt
N

aij
j =1

(|xi - xj |)+ fi (xi ), xi xi

i = 1 ...n.

(4)

This is our main equation and it fully defines the `coevolutionary' process of the system. In our model the evolution of a single species, by `climbing' to the local maximum of the fitness landscape, changes the landscapes of all those species which have an ecological link to it. Hence, the process described by (4) can be interpreted as a dynamical adaptation driven by ecological interactions. It results in a collective evolutionary motion of the species in the niche space. For the sake of simplicity in the following we will always assume that the evolution is neutral in the absence of all other species. This means that all external forces are set to zero, fi = 0. Moreover, we always assume periodic boundary conditions for the niche space.


226

The European Physical Journal Special Topics

Before we proceed with the analysis of the model some remarks are in order. First we note, that our formulation of a species' fitness bears some similarities to the seminal work by MacArthur [25, 26] who developed an elegant theory for the dynamics of n competiting species on a resource continuum. He assumed that each species i occupies a certain area of the resource spectrum according to an utilisation function ui (x), so that the competition strength between every pair of species is proportional to their niche overlap uij = ui (x)uj (x)dx. In this way the competition strength depends only on the species' distance on the resource continuum. Further, in the case of a Gaussian utilisation function ui (x) also the niche overlap uij is a Gaussian of the form Eq. (2) [25]. In a very similar way our model suggests that each species i = 1 ...n occupies a circular niche of radius in according to some (possibly Gaussian) utilisation function ui (x). The niche position xi = {xi,1 ,...,xi,m } is given by the centre of the utilisation function for each species and it is located at those conditions at which the species is most adapted to live. In this formulation it is natural to assume that the interaction strength between two species scales with the probability of interspecific interaction and we can thus interpret (xi , xj ) as the niche-overlap (xi - xj ) = ui (x)uj (x) dx. However, in contrast to the model by MacArthur our theory is not restricted to competition. Instead we assume that closeness in niche space will increase the strength of all types of species interactions, such as predation, mutualism, etc. While this approach is mathematically very convenient and allows us in a straightforward way to consider very different species in the same space , the plausibility of this basic assumption depends on the interpretation of the niche axis in the specific ecological system under consideration. For example, the position xi might quantify the habitat type where species i lives, the climate conditions it is best adapted to or the position along an environmental gradient. Typical examples include e.g., geographic altitude, humidity, light intensity, salinity or temperature. Clearly species, even if they are ecologically very different, will interact more strongly if they are similar in such an environmental space. In an estuary, for example, a two dimensional niche space could represent the salinity and turbidity of the water as environmental factors. Organisms living in the estuary might be best adapted to some positions along both salinity and turbidity gradients. Alternatively, the niche position might describe certain species traits, such as daytime of activity. Again, species with similar time of activity will have a greater chance to meet and interact, independent of whether the interaction is harmful or beneficial. Further, the niche space might also represent a chemical space of molecular locks and keys or signalling substances. Thus, we believe that our basic model assumptions should give an adequate description for a large class of environments and species traits. However, while this holds for some traits or environmental conditions, we are certainly aware that many ecological interactions cannot be described in this formalism. Most prominently the interaction strength will usually depend on many other factors than just the nichedistance. Also it is not clear that all the species in a community can be described in a single niche space. In a high-dimensional niche space one would assume that the interaction strength does not depend on the full distance but only on the distance regarding a certain sub-space or pro jection. Further it may well be possible that the functional form of (d) changes with the interaction type and for certain interactions does not decay with the distance d, and it may even depend on the spatial positions in niche space, (xi , xj ). Take, e.g. the body size (or related quantities such as skull length or beak size) which is a prototypical trait to determine the strength of interspecific competition. However, even though competition might be largest for two species of equal body size, this is usually very different for other interaction types. For example, in a predator-prey pair the predation rate should be maximal for a certain ratio of body sizes and also in parasite-host systems frequently strong interactions are achieved for drastic differences in body size of host and parasite. As another difference to the work by MacArthur and many others [25, 26] our model does not take species densities into account. This is an enormous simplification, because many evolutionary processes and parameters, such as mutation rates or extinction probabilities, depend on the population density. In this sense our model, at best, can only provide a sketch of the real evolutionary dynamics. On the other hand, this restriction makes model very general. While one certainly can generalise Eq. (4) to include population dynamics, e.g. of Lotka-Volterra type


Active Motion and Swarming
0.64

227
0.64

A y

0.64

B y

C

0.6

0.6

0.6

y
0.56

0.56

0.56

0.52

0.6

x

0.7

0.8

0.52

0.6

x

0.7

0.8

0.52

0.6

x

0.7

0.8

Fig. 1. Evolution of two species in a two-dimensional niche space for a pair of mutualists (A), competitors (B), or a predator-prey system (C). Solid and dashed lines show the tra jectory of the first and second species, respectively. The initial positions are labeled by open circles and the direction of motion is marked by arrows.

equations, such a procedure would include many additional assumptions (e.g. about analytic forms of saturation and density dependencies, relation of ecological vs. evolutionary time scales, etc.) and further give rise to technical problems that one would have to overcome, for example the avoidance of population explosions due to mutualistic interactions. All these problems and further choices do not arise in our model formulation, which makes our model very general and easy to analyse. Nevertheless we believe that many of the dynamic properties of our model will extend also to more involved formulations including population densities.

3 Two sp ecies mo del
In this section we consider the simplest case, namely the interaction of only two species, n = 2. In this case the matrix A of interspecific interactions has only two non-diagonal terms A= 0 a12 . a21 0 (5)

When only two species are present the gradient of the fitness function for each of the two species must point into the direction of the other species. Therefore, the resulting evolutionary dynamics is independent of the niche space dimension m and takes place on a one-dimensional straight line which connects the position of both species (see Fig. 1). Let z denote the position on the straight line connecting both species. Then the equations of motion are z1 = 1 a12 (|z2 - z1 |) z1 (6) (|z1 - z2 |), z2 = 2 a21 z2 where the fitness depends just on the niche distance d = z2 - z1 between the two species and might be given, for example, by our usual interaction function (d) = exp(-d2 /2 2 ). Suppose first, for simplicity that the mutation rates i = 1 and the interaction aij = ±1. Then the system Eq. (6) allows for three qualitatively different scenarios. In the first case of two mutualistic species the community matrix is given by the pair (a12 = +1,a21 = +1) and each species derives benefit from the other. In the evolutionary process both species "walk" on the fitness landscape (d) in search for a local fitness peak. Since (d) is supposed to be a monotonic decreasing function of |d|, for each of the two mutualistic species this peak is determined by the condition d = 0, i.e., it is placed exactly at the instantaneous position of their interacting partner. Hence the two species coevolve in such a way as to shorten their mutual distance in and they obtain the maximal fitness value when they encounter. At this point the system


228

The European Physical Journal Special Topics

reaches an equilibrium state. This is depicted in Fig. 1A which shows the tra jectories of two mutualistic species in a 2-dimensional niche space. We can interpret this dynamic outcome as a tendency of the mutualistic species to occupy similar habitats, i.e., they adapt to live in similar niches. The outcome of two-species competition (a12 = -1,a21 = -1) is shown in Fig. 1B. Here the interaction is prejudicial for both species, just in opposite to the previous case of mutualism. In the framework of our model competition means that the interaction is harmful if both species occupy similar niches, however it does not necessarily imply that they compete for similar resources. From Eq. (1) the local fitness landscape is such that each species has a fitness valley with the lowest value exactly at the actual position of its competitor. In the evolutionary process of climbing toward higher fitness both species move apart from each other and become increasingly separated in . The system reaches an equilibrium state when both species are located as far from each other as it is allowed by the boundaries. This behaviour of the two competiting species is well known from ecological literature as character displacement. In our case, using periodic boundary conditions, this means that the species are separated by a distance of half the system size. Finally, in a system of predator-prey or host-parasite type the interaction is given by an antisymmetric matrix with coefficients (a12 = +1,a21 = -1). In this case one of the species tries to avoid the contact, while its opponent benefits from the presence of its interacting partner and aims to shortens the mutual distance. Clearly the fitness function of the predator has a maximum exactly at the instantaneous position of its prey, while the prey fitness function has a trough at the position of its predator. Thus, the prey's favourable region is located as far from the predator as it is allowed by the system. In this way a directed movement in the niche space is generated and the pair moves uniformly with constant speed on a straight line in the direction given by the prolongation of their difference vector in . This dynamics mimics a coevolutionary arms-race (see Fig. 1C). Since the other species is also evolving, each of the two antagonists needs to "run for its life" in niche space just in order to keep its current fitness level (Red Queen Principle [7]). We want to stress however, that this arms race, where both species have precisely the same evolutionary speed, is only possible due to a fine tuning of parameters. In general, one would assume that this degeneracy is broken, either by one of the coefficients a12 or a21 being taken to be larger or due to some difference in the mutation rates i . In this case one of the species will obtain a larger velocity and the coevolutionary arms-race eventually comes to rest, either in the escape of the prey from or in its capture by the predator. To study this case we take the difference of both equations in Eq. (6). The niche distance d = z2 - z1 evolves as d = (1 a12 + 2 a21 ) (|d|). (7)

Usually we assume that the fitness peak is flat at the origin, (0) = 0. This means that the situation where both species have the same niche position z1 = z2 is always an evolutionary steady (but not necessarily stable) state. However, since (d) = 0 for positive distances d > 0 an arms race, where both species evolve with constant velocity d = 0, is only possible in the degenerate case, 1 a12 + 2 a21 = 0. Note that this can only happen for predator-prey or (+, -) type systems. In general, though, we expect that the system is degenerated and 1 a12 + 2 a21 = 0. Then, only two cases are left. If 1 a12 + 2 a21 > 0, we may speak of general mutualism. This situation might arise in two mutualistic species or in a predator-prey type system, if the prey evolves more slowly then the predator. Since we have assumed that the interaction strength is decaying with distance, (d) < 0, and has a local maximum at d = 0, the distance d(t) between the two species is always decreasing. This leads to a stop of the coevolutionary dynamics at full niche overlap. In the opposite case, of 1 a12 + 2 a21 < 0, we have a situation of general competition, which may arise in two competitive species or in a predator-prey system, if the prey evolves faster than the predator. Now, the distance is increasing in time, d > 0, and the two species move apart. From the ecological perspective this evolutionary outcome of a degenerated system, with the predator either catching the position of the prey or the prey escaping, possibly may lead to the extinction of prey or predator or both. However, in terms of our simple model this


Active Motion and Swarming

229

question cannot be answered since the model does not include abundances. In this context it is worth to note that here the predator fitness is always positive, whereas the prey's fitness is negative. This asymmetry stems from the fact that in the model other ecological relations and advantages of the prey, which do not depend on the niche position x1 (t) of its predator, are not taken into account even though they may contribute to the absolute fitness of the prey. In principle, we could make allowance for such factors by including constant external forces, fi (xi ) = ci , in the fitness function, Eq. (1). However, we refrain from doing so, because such terms will not change the resulting dynamics Eq. (4), which depend only the derivatives of the fitness landscape. This of course changes completely, if one aims to generalise the model in order to include species abundances, because then the absolute values of fitness become crucial. Here, it is just important to keep in mind that in our model absolute fitness levels are not well defined and it therefore does not make sense to weight the absolute fitness of different species against each other.

4 Anti-symmetric interaction matrix
As the analysis in the previous section has shown, in the absence of external forces fi (xi ) sustained motion in our model is only possible with degenerated predator-prey type interactions, while mutualisitic or competitive interactions favour stationary solutions. This corresponds to the believe of many biologists that predation is one main mechanism responsible for the evolutionary process [15, 16, 18, 20, 28]. In this study we are mostly concerned with coevolutionary motion in niche space. Therefore, for further analysis of our model with a larger number of species in the following we restrict us to the particular case of predator-prey interactions (a more in-depth investigation of the general case will be presented elsewhere). This restriction is not as severe, as it might seem at first glance. Predator-prey interactions are of extreme importance because they corresponds to food chains which are the basis of all ecosystems and, as we will show below, in our model they give rise to especially rich dynamics. Further, for sake of simplicity from now on we only consider the motion on a two-dimensional niche space, m = 2, with periodic boundary conditions. Mathematically, predator-prey type interactions relate to a system with an antisymmetric interaction matrix A, i.e., with elements aij = -aji . In the following we always assume that the mutation rates i are either identical for all species (or have been integrated already into the interaction matrix aij aij i ). One can easily show by taking the derivative of Eq. (4) that the antisymmetric system yields
n i=1

xi = 0. xi

(8)

This identity holds for all choices of the functional Thus the system is conservative and preserves phase Suppose first that there are only three species N general form 0a A = -a 0 -b -c and the equations of motion are given as

forms of the interaction strength (dij ). volume. = 3. In this case the matrix A takes the b c 0 (9)

dx1 =a (|x2 - x1 |)+ b (|x3 - x1 |) dt x1 x1 dx2 = -a (|x1 - x2 |)+ c (|x3 - x2 |) dt x2 x2 dx3 = -b (|x1 - x3 |) - c (|x2 - x3 |). dt x3 x3 (10)


230

The European Physical Journal Special Topics

Fig. 2. Basic three-species predator-prey configurations. (A) Three species chain, (B) feedforward loop, (C) feedback loop. Species are indicated as circles, the arrows correspond to the directed species interaction.

Further, the system preserves an integral of motion E = a(d12 ) - b(d13 )+ c(d23 ). Straightforward calculation shows that dE =0 dt (12) (11)

so that E (t) = const. This identity again can easily be verified by taking the derivative and using Eq. (4). The interaction matrix Eq. (9) allows for three different basic configurations (see Fig. 2). Following the classification of network motifs, introduced in [19] we will refer these as threespecies chain, feedforward loop and feedback loop. In the last case all species are equivalent, but in the first two configurations they are in a hierarchy. For brevity, we refer these species as top, middle and bottom species (or simply species 1, 2 and 3), respectively. 4.1 Three-sp ecies chain The three-species chain is the simplest configurations of a three species system (Fig. 2A). Given equal interaction strength the matrix of interspecific connections is given by the matrix (9) with a = c = 1 and b = 0 0 10 (13) A = -1 0 1 . 0 -1 0 This system might for example represent a three-trophic food-chain. In the evolutionary process the top species is chasing the middle species, which, in turn, tries to catch the position of the bottom species of the chain. Despite the similarity to the two species case, with N = 3 the dynamics is more rich and the behaviour of the system, depending on the initial conditions, can be completely different (see Fig. 3). Suppose, for example, a configuration in which the top and middle species are initially close together but are largely separated from the bottom species in niche space. If the niche overlap of the bottom species to the other two species is negligible then its fitness is not be affected by their position. Thus, the bottom species (3) nearly does not move in niche space, while the top and middle species form an almost isolated prey-predator pair (12), which moves in the fashion of an arms race with a constant speed, similar to the situation in Fig. 1C. Due to the periodic boundary conditions, this pair will eventually collide with the motionless bottom species. This is shown in Fig. 3A. At the time of encounter the system undergoes a switching of the middle species role from being persecuted by the top species to a chaser of the bottom species. The reason is that due to the attraction from the bottom species the middle species receives an additional acceleration and escapes from the top one. In this way a new arms race pair (23) appears and the former pursuer, the top species (1), loses the connection and stops evolving.


Active Motion and Swarming
1 0.8

231
0.5

A y

0.6

B

C
Switching

y

0.6 0.4 0.2

0.5

d23
0

0.4

Oscillations

0.3 0 0.2

x

0.6

1

0.2

0.3

x

0.4

0.5

-0.5 -0.5

d12

0

0.5

Fig. 3. Evolutionary dynamics in a three species chain. Depending on the initial condition the system exhibits either (A) switching in a collision of a predator-prey pair with an isolated species or (B) bounded, oscillatory motion. Solid line corresponds to the top species tra jectory, dashed line to middle and dotted line to bottom species. (C) Isolines of the first integral, Eq. (15), in the (d12 ,d23 )-plane.

Sooner or later, due the periodic boundaries, the new evolving pair hits the stopped top species and the reverse switching occurs in the similar manner, now, returning the middle species into the pair with the top one. Thus, the middle species turns successively to either persecuted or pursuer in different pairs. These collision and switching events resemble the scattering of elementary particles, which is well known from nuclear physics. In the three species case only two different scattering processes are possible (a pair of the top and bottom species (13) cannot persist because they have no direct link), which we symbolically can denote as 12 + 3 23 + 1 23 + 1 12 + 3. (14)

In contrast, if initially all species are located close to each other, then the system dynamics becomes very different. All species cannot leave some vicinity of each other and move together along oscillatory tra jectories around a central region in niche space (Fig. 3B). In this way beautiful quasiperiodic patterns of evolutionary movement are generated (see Fig. 4A). Further from Fig. 4B it becomes clear that the pairwise distances d12 (t) and d23 (t) are oscillating in anti-phase. This means that the middle species alternately is located close to the top or to the bottom species. It is therefore not possible to clearly distinguish between a bonded pair and an isolated species. But most importantly, in contrast to the two-chain in Fig. 1C, in a three chain the species do not evolve with equal speed. Instead, the coevolutionary motion of the three chain rather resembles a walking on three legs. The condition, which separates these two types of motion (switching or oscillation) can be found from the conserved quantity (11). In view of the interaction matrix (13) the integral of motion takes the form d2 d2 E = exp - 12 +exp - 23 . (15) 2 2 2 2 Apparently, each term in this sum is less than or equal to one. Hence, if E < 1 each of the exponential terms must be in the range [0,E ] and both distances are bounded from below, d > -2 2 ln E . In contrast, if E > 1 then each of the terms must be in the range [E - 1, 1] and therefore both distances are bounded from above, d < -2 2 ln(E - 1). Thus, the infinite motion, when one or another species leaves the group, is possible only if E < 1, while in the opposite case with E > 1 all species remain in a localised area. Fig. 3C shows the isolines, i.e., the lines with constant value of E , in the (d12 ,d23 )-plane. Since E is conserved, the system state must move along these one-dimensional lines. Switching orbits correspond to almost hyperbolic curves and oscillatory states correspond to the closed curves around the origin in the central area of the (d12 ,d23 )-plane. We note that this behaviour is reminiscent to the two-particle problem in a central force field in classical mechanics, where depending on the total energy of the system both particles form a bounded state or exhibit scattering.


232

The European Physical Journal Special Topics

Distances

Time

Fig. 4. Dynamic outcome in the three-species chain in the oscillating regime E > 1 (upper panel) and the switching regime E < 1 (lower panel). Plotted are typical tra jectories of the top species (solid line), the middle species (dashed lines) and bottom species (dotted lines) in niche space (left column) and the time evolution of the distances d12 (t) and d23 (t) (right column). In (C) only the tra jectories of the top and bottom species are shown to demonstrate that these two species do not enter a central circular region in niche space.

The fact that the motion takes place on one-dimensional curves at first glance would preclude the possibility of complex dynamics. However, as shown in Fig. 4C, this is obviously not the case. Clearly the evolutionary tra jectories follow complicated curves in niche space. This can be explained by the fact that the (d12 ,d23 )-plane does not represent the full phase space but is only a pro jection. Thus it is possible that the system exhibits complex and chaotic dynamics, even though it moves on a one-dimensional manifold in the (d12 ,d23 )-plane in Fig. 3. Clearly in Fig. 4C one can see the sharp edges in the tra jectories when species are scattered. Interestingly, the tra jectories of the top and bottom species never enter a central circular region in niche space. We do not have an explanation for this observation. The time series of the distances dij (t) in niche space show a very complicated dynamics (Fig. 4D) which suggests that the system is chaotic. This conjecture is confirmed from a calculation of the spectrum of Lyapunov exponents, i . In the three species chain we find that in the case of bounded dynamics E > 1 all 6 Lyapunov exponents are zero, i = 0. This confirms that in this case the dynamics is quasiperiodic. In contrast if E < 1 only 4 of 6 exponents are zero while the remaining two exponents are have finite value of opposite sign, 1 = -2 . This means that one Lyapunov exponent is larger than zero and the system is chaotic. The reason for this is obvious. The scattering as shown in Fig. 4A generates an instability and a sensible dependence on the initial conditions. Due to the periodic boundary condition the species are reinjected into the central area. As a consequence of this stretching (from scattering) and folding (from the periodic boundaries) chaotic dynamics can arise. Our numerical simulations in a system with an antisymmetric interaction matrix, A = -AT , show that the Lyapunov exponents always come in pairs with opposite signs, also for different network topologies and a larger number of species.


Active Motion and Swarming

233

4.2 Three sp ecies feedforward lo op We now turn to the next three-species configuration, the feedforward loop (Fig. 2C). Given equal interaction strength the interspecific connections are determined by the interaction matrix (9) with a = b = c = 1 0 11 A = -1 0 1 . (16) -1 -1 0 In comparison with the three-chain, this configuration has an additional link from the top to the bottom species. Thus, the top species is attracted by (and may bond to) both inferior species, whereas the bottom one avoids both superiors. The dynamic outcome again depends on the initial conditions. If the distance between the top and bottom species is large enough, then the extra connection can be neglected and the system demonstrates the same kind of evolutionary motion as the three species chain, namely either alternately switching of the middle species or quasiperiodic motion (see Figs. 5A and 5B). However, if initially the top and bottom species form a close predator-prey pair (13), then a new kind of dynamics appears. This is demonstrated in Fig. 5C which shows the species' tra jectories at the moment, when the top-bottom pair approaches the middle species. Now the bottom species experiences double pressure and becomes squeezed in from both sides by its superiors. As a result, in order to escape its opponents, it turns to the right or to the left relative to its tra jectory. The top species now experiences double attraction from both species and in principle may bond to either the middle or to the bottom species. However, since the bottom species is closer, the top species turns in the same direction and a switching of roles does not happen. Instead the top-bottom pair is repelled from the middle species and changes direction of motion. In the three species feedforward loop the integral (11) takes the form of E = exp - d2 d2 d2 12 +exp - 23 - exp - 13 . 2 2 2 2 2 2 (17) of motion. In lane. However, located one a approximately

So the extra connection yields a new term with negative sign in the integral contrast to Fig. 3C this should be depicted in a three-dimensional parameter p if we suppose a one-dimensional configuration, in which all three species are straight line in niche space, we can write d13 = d12 + d23 and the integral could be written in terms of d12 and d23 E = exp - d2 d2 (d12 + d23 )2 12 +exp - 23 - exp - . 2 2 2 2 2 2

(18)

In this form E depends only on two parameters and its isolines can be shown in a 2D plane (Fig. 5D). We found numerically that this approximation gives a good representation of the real dynamics. We can distinguish among the following dynamic regimes. First, the area of oscillatory motion, characterised by E > 1, with small pairwise distances in the central area in Fig. 5D. Note that this area is much smaller than the corresponding oscillatory area in the three species chain (Fig. 3C). Second, the area of switching in the range 0 < E < 1, with nearly hyperbolic isoclines similar to the corresponding picture in the three species chain. Finally, we find two diagonal areas which are characterised by negative values of the integral of motion, E < 0. In this regime the sum d12 + d23 = d13 , i.e. the distance between the top and bottom species, is much smaller than each value. This area corresponds to the repelling of the topbottom pair (13) by the middle species. Note, that the state with E 0 is located between the last two. If initially E 0, then the species move apart from each other and interact just slightly, even though they start quite close to each other. To summarise, the feedforward loop exhibits very similar dynamics as the three species chain. Large values of E lead to bounded states with quasiperiodic oscillations, whereas smaller positive values of E give rise to switching. We observe the same transitions as in Eq. (14). Thus even though, in principle, three predator-prey pair configurations are possible, switching just


234

The European Physical Journal Special Topics

A
0.6 0.7

B y
0.5

y
0.5

0.3 0.4 0.7 0.5

x

0.6

0.7

0.2

0.4

x

0.6

C
d

0.5

D

Oscillations Switching

0.5

y
0.3 0.1 0.2 0.4

23
0

Divergence
-0.5

Repelling
-0.5

x

0.6

0.8

d12

0

0.5

Fig. 5. Evolutionary dynamics in a three species feedforward loop. (A) Switching between pairs, (B) quasiperiodic oscillations around a central region, and (C) repelling of a (13) pair from the middle species. Solid line corresponds to the top species tra jectory, dashed line to middle and dotted line to bottom species. (D) Isolines of the first integral, Eq. (15), in the (d12 ,d23 )-plane.

leads to alternations between (12) and (23) pairs. The (13) however is associated with negative values of E and results in a new scattering process of repelling by the middle species. 4.3 Three-sp ecies feedback lo op The last antisymmetric three-species configuration is the feedback loop, shown in Fig. 2C. Assuming again equal interaction strength the interspecific connections are determined by the interaction matrix (9) with a = c = 1 and b = -1, giving rise to the matrix 0 1 -1 A = -1 0 1 . (19) 1 -1 0 In this configuration all species are equivalent, in that sense, that each of them can become a pursuer or a pursued species. This lead to similar patterns as in the three species chain, but with the main difference that now each of the species can change its role (see Figs. 6A and 6B). If the species are initially closely located, then as before they oscillate in a quasiperiodic fashion around one central point. Otherwise, if a predator-prey pair approaches an isolated species we observe switching. In contrast to the possible scattering reaction in the three-species chain in Eq. (14), the feedback loop gives rise to the following reactions 12 + 3 23 + 1 23 + 1 31 + 2 31 + 2 12 + 3. (20)

Interestingly, only three of the possible number of six scattering reactions are realised in this system. For example a 12 pair on collision with species 3 always forms a 23 pair and never a 31 pair.


Active Motion and Swarming

235
0.8 0.4

0.8

A

0.38 0.34

B

C

Switching

0.6

y
0.4

y

0.3

d23

0

0.26 0.22 0.4 0.6 0.8 1 0.25 0.3 0.35 0.4

-0.4

Oscillation
-0.8 -0.8 -0.4 0 0.4 0.8

0.2

x

x

d12

Fig. 6. Evolutionary dynamics in the three-species feedback loop. (A) Switching between pairs and (B) oscillations. Solid line is the top species tra jectory; dashed line is the middle and dotted line the bottom species. (D) Isolines of the first integral, Eq. (15), in the (d12 ,d23 )-plane.

In the three species feedback loop the integral (11) takes the following form E = exp - d2 d2 d2 12 +exp - 23 +exp - 13 . 2 2 2 2 2 2 (21)

Now the conserved quantity E can take values in the range from 0 to 3. Fig. 6C depicts the isolines of the feedback loop configuration, which correspond to a one dimensional configuration d13 = d12 + d23 similar to Eq. (18). In this figure the diagonal area corresponds to the middletop pair. But now this area is not separated from the other areas of infinite motion, as it was in the feedforward configuration. 4.4 N-sp ecies chain As an extension of three chain let us consider a n-species chain, i.e., a predator-prey hierarchy of n species, as is shown in Fig. 7. Assuming again symmetric interaction strength, the corresponding asymmetric interaction matrix is 0 1 0 0 ... 0 -1 0 1 0 ... 0 0 -1 0 1 ... 0 (22) A= . 0 0 -1 0 ... 0 ... ... ... ... ... 1 0 0 0 0 -1 0 It is easy to verify, that the integral of motion E in this case takes the following form
n-1

E=
i=1

exp -

d2 i+1 i . 2 2

(23)

Further we find that the number of zero Lyapunov exponents is always four, which corresponds to the symmetries Eq. (8) and Eq. (11), while the remaining 2n - 4 exponents come in pairs with opposite sign. The resulting evolutionary dynamics in niche space are shown in (Fig. 8). In dependence of the value of the conserved quantity E in Eq. (23), we can distinguish between n - 1 different dynamic states of the system. If E < 1, then only one predator-prey pair can persist at one time instance and we observe switching. If k < E < k + 1, where k is an integer number, then the model always contains k predator-prey pairs, which exhibit complicated patterns of switching (Fig. 8A). And finally, a large value of E > n - 1 means that all distances are sufficiently small and all species moves in an intricate collective motion, which resembles a swarming pattern (Figs. 8B and 8C).


236

The European Physical Journal Special Topics

Fig. 7. Schematic representation of a n-chain.

We have found that the dynamics of such a swarm-like, dense collection of species depends on the parity of the species number n. For example, a single species in the absence of others (n = 1) does not evolve and always remains at the same position in niche space, whereas a prey-predator pair (n = 2) moves uniformly and linearly, as shown in Fig. 1C. A chain of three species (n = 3) again, if initially located close enough to each other, does not show a mean net displacement and just rotates around a central area (Fig. 5), whereas a chain of 4-species (n = 4) on average moves as a whole along one direction (not shown). We conjecture that this alternating behaviour of the collective movement with the species number continues for an arbitrary number of species of the chain (always given the proper initial positions). If a chain consists of an even number of species (e.g, n = 2, 4, 6,...) its centre of mass shows a directed movement and the chain as a whole moves in a linear and uniform way through , in a way that resembles the motion of a two species predator-prey pair. In contrast all chains with an odd number of species (e.g, n = 1, 3, 5, 7,...) just oscillate around one central location and statistically stay at the same region in . While we do not have an hard evidence of this observation, we have verified numerically that this dependence of the collective motion on the parity of the chain is true for chains up to 20 species. In the other case, for lower values of E if the species are not close enough initially, then the minimum niche overlap among them is not achieved and the chain will not evolve in as a whole. Instead it decays into evolving sub-chains conformed by species that are close enough. These sub-chains follow the same rules of motion as noted previously. For example, the sub-chains may collide, with the usual result of species-switching ongoing with changes in the direction of motion, much like it was discussed in detail in the 3 species example. Unfortunately we cannot formulate a law (an analog of conservation of momentum), which would show the direction of motion after collision. This does not hold only for linear n-chain topologies, but for arbitrary antisymmetric interaction matrices. Further (again without formal proof ) we have observed the following rule for species switching, when a fast evolving n-chain collides with an isolated species. Obviously, species-switching can only happen if the initially isolated species has an interacting link with the oncoming chain. However, we have found that this link must be such that it is lower in level than the lowest species in the chain or higher in level than the top species in the chain. In other words, the single interacting species must be lower or higher in hierarchy than the lowest or highest species of in the coevolving chain, respectively. If this is the case, then the swapping always takes place. For example, if an arms-racing pair meets an isolated species that
0.62

A
0.6

1.4

B

c

0.58 1

y

0.4

y
0.6

y0.54
0.5

0.2

0 0.2 0.4

0.2

x

0.6

0.8

1

0.2

0.6

x

1

1.4

0.2

0.25

0.3

x

Fig. 8. Evolutionary dynamics in a linear n-chain. (A) Dynamics is separated into several sub-clusters which exhibit switching of pairs. (B) Swarming motion with a directed collective motion in a chain of 10 close species. (C) Swarming motion without net displacement in a chain of 9 close species.


Active Motion and Swarming

237

can be positioned in the lowest level of the chain, then the top species is out-raced and the middle species bonds with the former isolated species by forming a predator-prey pair in which the middle species is chasing the incoming new species. However, the evolutionary dynamics is entirely different if the pair encounters a species that can be placed in the middle of the chain hierarchy (the feedforward loop). In this case the switching does not occur and, on the contrary, the incoming pair is repelled by the isolated species.

5 Discussion
In this article we have presented a simple conceptional model for the coevolutionary dynamics of ecological interacting species. Similar to other recent studies [24, 27] we propose that the species evolve by changing their positions relative to the other species in a common niche space, which might for example represent a hypercube of certain phenotypic traits or environmental conditions to which the species are best adapted. Depending on the interaction type (e.g., mutualistic, competitive, ...) each species can have either positive or negative impact on the fitness of its interacting partners. Our basic assumption is that the strength of the interaction between two species, independently of the interaction type, will be determined by the species' niche overlap, i.e. their pairwise distance in niche space. Thus, the main idea of our model is to separate the contributions to the evolutionary fitness of a species into the product of a time-independent community matrix, which determines the signs and the type of the species interaction, and time-dependent distance functions which are responsible for the strength of the interaction. Using numerical simulations we have investigated the resulting evolutionary motion in a twodimensional niche space. We have explored in detail all antisymmetric three-species network configuration as well as a range of linear n-chains, from a single pair up to 20 species. These investigations revealed that the evolutionary dynamics increases in complexity and richness with the number of species in the model. In the simplest case of only two interacting species we observe a rather straightforward model outcome, where the motion takes place on a straight line. If the community matrix is symmetric the system exhibits either a convergence/divergence equilibrium, where two competitors move in opposite directions to avoid the harmful influence of the other (character displacement), while two mutualists approach until they encounter in one point. In contrast, in a system with an antisymmetric community matrix, e.g. a prey-predator pair, both species show a linear and uniform movement in the same direction (coevolutionary arms race). More complicated dynamics appear if more species are present in the community. Here we were interested in a system with an antisymmetric community matrix, where we have shown that the flow of the species in phase space is conservative and preserves a conserved quantity. Nevertheless the evolutionary dynamics is characterised by seemingly turbulent dynamics, where the continuous adaptive changes in the evolutionary landscape drive the species into an intricate collective swarming pattern. Despite this complexity, we were able to identify several dynamical rules. First, we found that as a prerequisite for coevolution the species must be sufficiently close in niche space, which means that a minimal niche overlap must exist. As a consequence, depending on the initial conditions the system separates into an ensemble of quasiindependent subsets of species which form more or less linear chains. Depending on the parity of their species number these sub-chains either oscillate quasiperiodically in a localised domain or evolve with near constant velocity in a preferred direction. Although it is widely accepted that the evolutionary process in biology in general is not directed (for a profuse discussion see [17]), the unexpected formation of such moving species-swarms in our model shows the emergence of directed evolution. That is, the `evolutionary direction' of the collective system arises as an emergent property of the individual rules (e.g., uphill fitness climbing) and the topology of the interactions. We note, though, that the evolutionary path of individual species in such a moving chain is not uniform but is characterised by alternating or antiphase oscillations of next-neighbour distances in the chain. Thus, our model predicts that the directed evolutionary motion of n species does not resemble the uniform movement of a n-particle rigid body but would be better described as `evolutionary walking' similar to a n-legged millipede.


238

The European Physical Journal Special Topics

Furthermore, when such sub-chains collide in their evolutionary motion, we observe chaotic scattering or repelling with drastic changes in the direction of movement, ongoing with an exchange (or switching) of individual species between sub-chains, reminiscent to the scattering of elementary particles. Further, for a small number of species the evolutionary dynamics exhibits a staring analogy to the central force problem in classical mechanics, where depending on the total energy of the system the particles form a bounded state or exhibit scattering. We want to stress that these evolutionary collision and scattering events are a non-trivial result of our model and result in punctuated, drastic changes for an ongoing evolutionary arms race within a subchain, which may either abruptly change its evolutionary direction, vanish completely or form a new different race with selected species from the collision partner. As a consequence any species in a coevolutionary arms race potentially can be `out-raced' by the influence of other interacting species. Even when we are aware that the presented model is very simple and holds strong assumptions, not least the fact that we do not include population densities, we believe it well worth to be considered, given the non-trivial features it produces to the elaboration of more realistic biological models and theory.

We thank A. Pikovsky and M. Zaks for helpful discussions. This work was supported by German Volkswagen Stiftung, DFG (SFB 555) and the Venezuelan foundation Fundayacucho.

References
1. J.M. Chase, M.A. Leibold, Ecological Niches: Linking Classical and Contemporary Approaches (University of Chicago Press, Chicago, 2003) 2. S.P. Hubble, The Unified Theory of Biodiversity and Biogeography (Princeton University Press, Princeton, 2001) 3. J. Grinnell, Auk 34, 427 (1917) 4. C. Elton, Animal Ecology (Sidgwick and Jackson, London, UK, 1927) 5. R.S. Miller, Adv. Ecol. Res. 4, 1 (1967) 6. G.E. Hutchinson, A Treatise on Limnology (Geography, Physics and Chemistry, 1957) 7. L. Van Valen, Evol. Theory 1, 1 (1973) 8. N.C. Stenseth, J. Maynard Smith, Evolution 38, 870 (1984) 9. D.A. Rand, H.B. Wilson, Proc. R. Soc. Lond. B 253, 137 (1993) 10. M.F. Dybdahl, C.M. Lively, Proc. R. Soc. Lond. B 260, 99 (1995) 11. D. Ebert, W.D. Hamilton, Trend. Ecol. Evol. 11, 79 (1996) 12. E. Szathm´ , I. Scheuring, C.S. Hegedus, G. N´ ary ~ emeth, I. Moln´ G. Vida, Coenoses 5, 130 (1990) ar, 13. R.V. Sol´ R. Ferrer, I. Goz´ e, alez-Garc´ J. Quer, E. Domingo, J. Theor. Biol. 198, 47 (1999) ia, 14. F. Dercole, R. Ferrier` A. Gragnani, S. Rinaldi, Proc. R. Soc. B 273, 983 (2006) ere, 15. R. Dawkins, The Extended Phenotype (Oxford Univ. Press, Oxford, 1982) 16. S.J. Gould, Ever Since Darwin: Reflections in Natural History (Norton, New York, 1977) 17. S.J. Gould, The Structure of Evolutionary Theory (The Belknap Press of Harvard University, Cambridge, 2002) 18. C.K. Levy, Evolutionary Wars (W.H. Freeman, New York, 1999) 19. R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, U. Alon, Science 298, 824 (2002) 20. G.J. Vermeij, Annu. Rev. Ecol. Syst. 25, 219 (1994) 21. C. Hauert, S. De Monte, J. Hofbauer, K. Sigmund, Science 296, 1129 (2002) 22. E. Kisdi, F.J.A. Jacobs, S. Geritz, Selection 2, 161 (2001) 23. U. Dieckmann, P. Marrow, R.J. Law, J. Theor. Biol. 176, 91 (1996) 24. H.C. Ito, T. Ikegami, J. Theor. Biol. 238, 1 (2006) 25. R. MacArthur, R. Levins, Am. Nat. 101, 377 (1967) 26. R.M. May, R. MacArthur, Proc. Natl. Acad. Sci. 69, 1109 (1972) 27. A. Branstrom, U. Diekmann (unpublished) 28. P.A. Abrams, Annu. Rev. Ecol. Syst. 31, 79 (2000)