- Email: [email protected]
The dynamics of group formation
ELSEVIER
The Dynamics of Group Formation SHAY GUERON Department of Mathematics, Technion-- lsrael Institute of Technology, Haifa, 32000 Israel AND
SIMON A. LEVIN Department of Ecology and Evolutionary Biology, Princeton University, Princeton, New Jersey 08544-1003 Received 6 January 1994; revised 10 April 1994
ABSTRACT A general continuous model is presented for animal group size distribution. Attention is restricted to a fixed size population divided into groups of various dynamic sizes, but the approach extends easily to populations of variable size. The basic idea is to relate group size distribution to two functions, the (densitydependent) rates of fusion and fission. These functions can be estimated from data and can ultimately be related to the behavior of individuals and the dynamics of groups. For various functional forms, the stationary distributions of group sizes are sought. In several prototype cases, the stationary distribution has a peak value, the "most frequent group size," which emerges endogenously from the dynamics. The authors determine when such a peak emerges and more generally show the existence and uniqueness of the stationary distribution. Stability of stationary solutions is discussed. Progress is shown, but a general treatment remains refractory.
1.
INTRODUCTION
The problem of animal grouping is one of the most fundamental ones in biology. The nonuniform aggregation patterns of organisms ranging from small invertebrates (see, e.g., [6, 12]) to large vertebrates (e.g., [14, 16, 17, 22]) have both ecological and evolutionary significance, and the tendency to aggregate is under strong evolutionary control. The causes and consequences of grouping represent fascinating subjects for theory and will be discussed elsewhere (see, e.g., [8, 15]). In this paper, we restrict attention to relating observed grouping patterns to the processes of aggregation (fusion) and splitting (fission). In theory, both of these processes can be estimated from data and ultimately can be related to individual behavior. For example (see [11]), the decision of MATHEMATICAL BIOSCIENCES 128:243-264 (1995) © Elsevier Science Inc., 1995 0025-5564/95/$9.50 655 Avenue of the Americas, New York, NY 10010 SSDI 0025-5564(94)00074-A
244
SHAY GUERON AND SIMON A. LEVIN
an individual to join or not join a group of size n, as well as the decision by individuals in a group of size n to permit or not permit an additional individual to join can be represented in the fusion rate of groups of size 1 and groups of size n. Even the fundamental problem of pair formation/disassociation is a special case of a fission/fusion process, albeit one for which a two-sex model would be more appropriate. Some models for animal aggregation suggest that the underlying grouping mechanism involves individuals "seeking target density" (STD). In such models the individual is assumed to prefer a target number of neighbors within its perception/detection range and responds to its perceived environment accordingly. For example, based on STD, GriJnbaum [7] derives an integrodifferential equation for swarming. In a complementary approach, Gueron and Liron [9] discuss a situation where movement is influenced by a combination of chemotactic attraction and social behavior. They show that STD does not yield grouping, at least not in the restrictive sense of traveling wave solutions, whereas other social behaviors can. STD and other models that rely on local cues do not deal with group size distributions explicitly. A central question is how group size distributions in nature arise from the rules individuals follow and the rates for fusion and fission. Most important, the modal group size may be an epiphenomenon, not one preferred or even perceived by individuals. The approach presented here is not entirely novel, either within the framework of population biology or elsewhere. Cohen [3] has developed a fundamental treatment of the group size distribution problem but restricts transformations of groups to those that cause the group to increase or decrease by a single individual. Other approaches include that of Anderson [2], who presents and analyzes a linear model for the dynamics of change. Finally, a sophisticated literature exists within the study of polymers (e.g., [1, 4, 5, 10, 18-21]), in which models very similar to those we analyze herein are discussed. 2.
DEVELOPMENT OF THE BASIC MODEL
We regard the continuous model presented here as a limit of a discrete model (see, e.g., [11]). The continuous treatment facilitates calculation and introduces a number of delicate points of interpretation. For example, group sizes are allowed to vary over [0,~) (but the distribution is constrained to have an exponentially decaying "tail") although the population is of a fixed finite size. In the remainder of this section we eschew delicacy and restrict attention to the continuous model. We consider a population of total size P, distributed into groups of size x (x~[0,~)), where f(x,t) denotes the density distribution of
DYNAMICS OF G R O U P F O R M A T I O N
245
groups of size x at time t. For convenience, we drop the explicit notation for the time dependence. Instead, we distinguish time dependence by a tilde [e.g., f(x) stands for f(x, t)]. Assuming that all of the moments (about 0) of f(x) exist, we denote
Mn = fo z"f( z)dz =
the nth moment about O.
We introduce a special notation for the first two moments, since these have important biological interpretation: ~c
(~ P
=
M0
= fo f( z)dz =
=/~1 =fo
the number of groups,
zf( z)dz = the total population size.
Recall that P is fixed in our model, although it would be an easy extension to relax this assumption (see [11]). In the following sections we discuss the stationary distributions [denoted f(x)] and their associated moments (Mn). These are, by definition, time-independent and hence the tilde is omitted. The model rests on two density-dependent functions, qJ and ~b, describing rates of fusion and fission of groups. In particular,
tO(x,y) =
the probability per unit time for two groups of sizes x and y to merge. ~h(x, y ) = the probability per time for a group of size x to fragment and generate groups of size y and x - y (we ignore breakage into multiple fragments). We also introduce
p(x) =
the probability per unit time for a group of size x to break into two fragments.
Note that p(x) is not independent of the other functions but is simply (because of double counting) half the integral of 4, over all possible group sizes smaller than x [see (2.4) below]. Clearly, we can also introduce a function Q to satisfy
qb(x,y) =p(x)Q(x,y),
(2.1)
where Q(x, y) is the probability that, given that a group of size x breaks into two fragments, one of them is of size y (and the other is of size
x-y).
246
SHAY GUERON AND SIMON A. LEVIN
The quantities fb(x,y), qJ(x,y), p(x), and Q(x,y) must satisfy the following requirements:
(2.2)
(by definition),
x
(2.3)
(fragments are smaller than the "mother" group),
p(x) =3
(2.4)
(by definition), so that
foQ(X, z)dz = 2
(2.5)
qJ( x,z) = qJ( z,x)
(2.6)
[provided p(x) 4: 0] and
(symmetry). 3. A GENERAL INTEGRODIFFERENTIAL BALANCE EQUATION AND ITS STATIONARY DISTRIBUTION 3.1. TIME-DEPENDENT INTEGRODIFFERENTIAL BALANCE EQUATION The time evolution of the density distribution of groups is the result of the balance between fission (fragmentation of groups into smaller groups) and fusion (merger of groups into larger groups). It can be described by the following balance equation (see also [1]):
Ot
p(x)f(x)- fof(X)f(z)
1£.
+ -~
f(y)f(x-y)O(y,x-y)dy+ III
jxf ( y ) 4 ) ( y , x ) d y . IV
(3.1)
The four terms on the right-hand side of Equation (3.1) account for (I) fragmentation of groups of size x, (II) amalgamation of groups of size x
247
DYNAMICS OF GROUP FORMATION
with other groups, (III) generation of groups of size x by merging of two smaller groups; and (IV) fragmentation of larger groups into fragments, one of which is of size x. 3.2. THE EQUATION FOR THE STATIONARY DISTR1BUTION
A stationary distribution is a solution of Equation (3.1) that satisfies
Of(x) 3t
O.
(3.2)
The equation for this stationary distribution, f(x) (tilde omitted), is
- p(x)f(x) +-il
- f ( x ) foxy(z)
z) dz
f; f ( y ) f ( x - y ) O ( y , x - y ) d y +
ix~f ( y ) 4 , ( y , x ) d y = O .
(3.3)
We investigate the existence, uniqueness, and stability of such a stationary distribution. Furthermore, to allow comparison with data, we solve for this distribution in special cases. 4.
A MODEL TEMPLATE AND T H R E E TYPICAL EXAMPLES
For the remainder of this paper, we assume that Q(x, y) is independent of y; that is, the conditional distribution for fragments is uniform: 2
Q(x,y) =x'
(4.1)
which implies 2
ch( x,y) = x P ( X ) .
(4.2)
This somewhat restrictive assumption represents the most reasonable null case, but it will be important to consider other cases in future work. Our model is restricted to cases where the rates at which groups undergo fusion and fission depend only on the groups' sizes. Therefore, the (time-independent) functions O(x,y) and p(x) reflect the "social dynamics" of the population. We further assume that the symmetric function O(x, y) is separable, so that
qJ( x, y ) = aa( x )a( y ), where a(x) is some nonnegative function.
(4.3)
248
SHAY GUERON AND SIMON A. LEVIN
The resulting equations can be investigated numerically, even though they are resistant to analytical investigation. It is instructive, however, to consider the special case
p(x) = ~Sxa(x),
(4.4)
for which an explicit solution is feasible. As shown in the next section, this case carries substantial biological interest. Together with Equation (4.2), Equation (4.4) implies
~b(x,y) = 2/3a(x).
(4.5)
With this, the time evolution equation (3.1) becomes 0f(x,t)
Ot
Bxa( x)f( x ) - aa( x)f( x) rosa( z)f( z)dz + ½ a ~ a ( y ) f ( y ) a ( x - y ) f ( x - y)dy +2#~a(y)f(y)dy,
(4.6)
whereas the equation for the stationary distribution (3.3) translates to
- ~ x a ( x ) f ( x ) - aa(x)f(X)fo~a(z)f(z)dz + ½afo a ( y ) f ( y ) a ( x - y ) f ( x - y)dy + 2~
Jo"a ( z ) f ( z ) d z = O.
(4.7)
At this point, any specific case of the model is defined by the choice of the function a(x) and the parameters a and /3. 4.1. THREE 7)'P1CAL EXAMPLES We focus our attention on three typical examples, where the function
a(x) is decreasing, increasing, or constant. As shown below, these yield three different types of stationary distributions. Of course, other variants may easily be adopted.
Case 1 a(x) = 1.
(4.8)
In this case group fusion is independent of group size; that is, the probability for a group of size x to merge with a group of size y is
D Y N A M I C S OF G R O U P F O R M A T I O N
249
constant (independent of x and y). Furthermore, group fission [represented by p(x)] is proportional to group size, again quite a reasonable null hypothesis. Consequently, ~0(x,y)= ~, p ( x ) = / 3 x , and ~b(x,y)= 2/3 (c~ and /3 are positive constants). Case 2
(4.9)
a ( x ) = x,
which implies that ~(x, y) = axy, (o(x, y ) = 2/3x, and p ( x ) = fix 2, This accounts for grouping behavior in which large groups are more "attractive" or likely to collide and merge with others than are smaller ones, in proportion to their size. Furthermore, larger and larger groups become increasingly unstable, so that the fragmentation rate increases nonlinearly with group size. Case 3
(4.10)
a(x) = 1/x,
which implies that ~b(x,y)= a / x y , ¢ ( x , y ) = 2 / 3 / x , and p ( x ) = / 3 . Thus, groups break at a rate independent of their size, for example, due to external disturbances that affect all groups uniformly. Furthermore, the functional form of ~O means that solitary individuals, or individuals in small groups, do not remain in such groups for long. Possibly, due to the fact that it is disadvantageous to be in small groups, small groups merge at high rates to form larger ones. 5.
IDENTITIES, R E C U R R E N C E RELATIONS, AND CONVERGENCE PROPERTIES
5.1. IDENTITIES Here we prove two identities that are used later. Let f ( x ) be a function, defined for x~[0,~), that decays to 0 exponentially as x--* ~. Then all of its moments about zero (denoted M,) exist, and we obtain
=L
I(y)dyf,
+ yl"f(z)dz n
=
C.MkM,,_ k k=0
(5.1/
250
SHAY GUERON
A N D S I M O N A. L E V I N
and
fo xn d~ fx f ( z ) dz = - n +1l 5.2.
( z ) dz Sozx" d~ fo~z " + l f ( z ) d z = n +1l M,+ 1. (5.2)
DERIVING R E C U R R E N C E R E L A T I O N S FOR THE M O M E N T S
We define the auxiliary function
(5.3)
e(*) = a ( x ) h x ) and denote its moments about zero by Nn = f0~xn~(x) dx
(5.4)
(assuming that N,, M,, N,, and M n exist for all natural n and for all t). A recurrence relation between the moments of f ( x ) and g(x) can be obtained by multiplying both sides of Equation (4.6) by x ~ and integrating from 0 to ~ with respect to x. Applying the identities (5.1) and (5.2), we obtain
d& .
dt
.
- 1., + lf. l ~n N
aNolqn +½a ~
C nkN" k N "n _ k
k=O n-1
=-
n - 1 N# - + i +½a E /3 ~---~
k" "
CnNkNn_ k .
(5.5)
k=l
For n = 0, one must use only the "explicit" (first) equality in Equation (5.5). For n = 1, Equation (5.5) is an identity with both sides equal to zero. This is a manifestation of the conservation of the total population size (/~a). Note that for our three special cases defined above, we have /Vn = h4n
[for a(x) = 1],
(5.6a)
N~ =/~n÷l
[for a(x) = x],
(5.6b)
[for a(x) = 1/x].
(5.6c)
N~÷1 = / ~ 5.3.
C O N V E R G E N C E P R O P E R T Y FOR THE CASE a ( x ) = 1
The case a(x)= 1 is the only case where the system of coupled differential equations (5.5) yields a straightforwardly solvable equation, at least with our skills. We derive an equation for the evolution (in
DYNAMICS OF GROUP FORMATION
251
time) of the number of groups, (~ ( = M0), by using the identity (5.5) for n -- 0. From (5.6) we have M n -- ]Qn, which gives dG
1
dt = - T A G
- 2
+tiP
(P=2Q,).
(5.7)
The solution of Equation (5.7) can be found analytically by substitution of
U
=½.d.
(5.8)
This yields
mexp{[(1/2)tap]'/zt}_nexp{_[(1/2)~ap]'/2t} (~(t) = i -2t-Pa m e x p { [ ( 1 / 2 ) f l a P ] '/2 t} + n exp{- [ ( 1 / 2 ) f l a P ] 1/2t} (5.9) where the ratio m / n is determined by the initial conditions. It is now easy to see, as was originally shown by Aizenman and Bak [1], that regardless of initial conditions the number of groups converges to the limit lim (~(t) = i 2fl--e . 6.
(5.10)
EXISTENCE AND UNIQUENESS OF THE STATIONARY DISTRIBUTION
Although we cannot solve the more general problem or prove stability, we can establish the existence and uniqueness of the stationary distribution for Equations (4.1)-(4.7). Note that these are more general than the three special cases mentioned earlier. To that end, it is important to notice first that if one defines the auxiliary function g(x) by
g(x) - a ( x ) f ( x )
(6.1)
and writes Equation (4.7) in terms of g(x), it is reduced to the form of Equation (4.7) written in terms of f(x) for the case a(x) -- 1. Therefore, we conclude that showing the existence and uniqueness of the stationary distribution for the special case a(x) = 1 implies the same result for the more general case.
252
SHAY GUERON AND SIMON A. LEVIN
We denote the Laplace transform of f(x) by
fo e-SZf( z)dz - F( s).
(6.2)
F ( 0 ) = G.
(6.3)
zc
Sa( f( x) ) = For s = 0 we have
To solve the equation for the stationary distribution we take the Laplace transform of Equation (4.7) and obtain
/3F'(s)- o~GF(s)+½o~FZ(s)+2/3G - 2 / 3 F ( s ) = 0 . s
(6.4)
s
Equation (6.4) is a first-order O D E with a singular point at s = 0. For every s o > 0, this equation is regular on the interval [s0,~). As such, all of its solutions on this interval belong to a one-parameter family, namely, /3 1 F ( s ) = 2 S (s-~--d).
(6.5)
Continuing a solution from the above family to s ~ [0, ~) determines the parameter d uniquely, and it is easy to verify that the only solution of Equation (6.4) on the interval s ~ [0, ~) is s +2(/3/o~G)
"
(6.6)
Transforming F(s) back to the real domain, we obtain the stationary distribution, which is uniquely determined by G:
f( x) = 2 ~ e -2(tJ/'~c)x OL
(6.7)
Now, instead of fixing G, one can equivalently choose P as the normalization factor. For a given value of P, the stationary solution is
and the number of groups (G) is related to P by
G = ~/-~~
(6.9)
DYNAMICS OF GROUP FORMATION
253
Returning to the general case we write the stationary solution e
where the value of h is determined by the equation P=2--~
a - z~
e AZdz.
(6.11)
This result is valid for any positive a(x) with algebraic asymptotic behavior at x ~ ~ (to ensure that all moments exist). Singularity occurs when a(0)= 0 [e.g., for a(x)= x], although the integral may still converge. Consequently, to ensure that the population (P) remains finite, we must restrict the choice of a(x) to satisfy lim0 -a(7 -x ) > 0. ,-~
(6.12)
All three special cases satisfy these conditions. Figure 1 displays the shapes of the relationships (6.10) for the three special cases, with a, /3, and P all set equal to 1. 7.
I N T E R P R E T A T I O N OF T H E STATIONARY SOLUTIONS
From Equation (6.10) it is easy to see that the properties of the resulting stationary solution depend only on a(x). Also, it is easy to integrate the expressions for f(x) and xf(x) over any desired interval [xl,x 2] to obtain the number of groups whose sizes are in this interval and the associated total population. Equation (6.10) demonstrates that for every choice of a nondecreasing function a(x) [e.g., a(x)= 1 and a(x)= x] the resulting stationary distribution is a monotonically decreasing function of x. For a(x)= 1 we obtain an exponential distribution of group sizes in which the most frequent groups are the "small ones." The same conclusion is true for a(x)= x, where the resulting distribution f(x) is singular at x = 0. In this case we have small groups piling up near x = 0 and larger groups becoming less abundant. On the other hand, we conclude that any choice of a decreasing function a(x) [e.g., a(x) = 1/x] yields a group size stationary distribution that has at least one local maximum. This maximum is the modal group size. It follows, not surprisingly, that a population can display a modal group size without having this quantity built explicitly into the dynamics.
254
SHAY GUERON AND SIMON A. LEVIN
2.5-
0
0.5
1
1.5
2
2.5
3
3.5
(a~ FIG. 1. The shapes of the stationary solutions in (6.10) for the three special cases (a) a ( x ) = 1, (b) a ( x ) = x, (c) a(x)= 1/x. In all three plots, a, /3, and P are set equal to 1. The resulting values of h are (a) h. = Vr2, (b) h. = 2, (c) h = 41/3.
DYNAMICS OF GROUP FORMATION
ol 0
I 0.5
255
t
t
1
1.5
x
(b) FIG. 1. Continued.
256
SHAY GUERON
A N D S I M O N A. L E V I N
iT
0.8-
0.6-
0.4-
,
\ 0.2-
\
\
\\ \
\
h
2 x
(c) FIG. 1. Continued.
L
3
t
4
DYNAMICS OF GROUP FORMATION 8.
257
STABILITY ANALYSIS
In this section we present work in progress; the full problem remains impenetrable. We investigate the local stability of the stationary solutions for the cases a(x) = 1 and a(x) = x. As demonstrated below, this leads to a nontrivial eigenvalue problem, for which we have only a partial solution. Our hope is that others will be inspired to complete the proof, and hence we present these results as open questions. 8.1. EIGENVALUE PROBLEM FOR A GENERAL PERTURBATION We start with the dynamical balance PDE (3.1) written in terms of ff,(x) - f(x)a(x). With no loss of generality we assume ~ =/3 = 1 (this can be achieved by time and group size rescaling) and obtain
~,(x) a( x)
x~,( x)
'S/
+g We consider 2e-AX/a(x),
-
~,(x)f° g(y)dy
g(x-y)g(y)dy+2
a perturbation
f(x) =
on the
/;
g(y)dy.
stationary
2e :ax + h ( x ) a(x) '
solution
(8.1) f(x)=
(8.2)
for some (time-dependent) function h(x). Conservation of the population size restricts the initial perturbation /~(x)t=0, and the integral constraint
t
X~a(x)],=o
=o
(8.3)
must hold. To allow for linearization, we start with /Tt(x),=0 < O.
(8.4)
258
SHAY GUERON AND SIMON A. LEVIN
We linearize Equation (8.3) by omitting the two quadratic terms
fo~(x-y)~(y)dy,
[~(X)fo~h(y)dy
(8.5)
and obtain the linear equation, which, after reorganization, can be written as
-
+ 2e-AX~t(y)e ~' dy + 2fx~[t(y)dy.
(8.6)
Assuming that/~ can be expanded in terms of eigenfunctions, we use the standard separation of variables
tl = emtk(x)
(8.7)
and obtain an equation for k(x):
2 m k(x) a(x) =_ xk(x)--xk
_2e_.Xfok(y)dy
+ 2e-~Xfok(y)eXY dy + 2fx k(y)dy.
(8.8)
Defining the linear operator 5r by k ( x ) ) - - x k ( x ) - -2 2k
_2e_~fo~k(y)dy ~z
+2e-~Xfok(y)eXYdy+2fx k(y)dy,
(8.9)
we regard the problem of solving Equation (8.8) as being one of finding "weighted" eigenvalues of the operator 5r,
~( k( x) ) = m k( x) a(x) [where the weight function is complex.
1/a(x)].
(8.10)
Note that m is in general
DYNAMICS OF GROUP FORMATION
259
We point out that 5r is a "conserving operator." Multiplying both sides of Equation (8.8) by x and using the identities (5.1) and (5.2), we have ~c
fo x~(k(x))dx = 0.
(8.11)
As a result, for any eigenfunction k(x) that satisfies (8.10), the perturbation component of the form (8.7) conserves the population size automatically. The integral equation (8.8) can be converted to a second-order ODE. This is achieved by differentiating it twice with respect to x (assuming the necessary smoothness), which yields
m[a--~
+
+ Am k(x) a(x) ] + 3 A k ( x )
(Ax+4)k'(x) 0,
(8.12)
where a prime denotes differentiation with respect to x. For a given choice for the function a(x), by expanding the solution as a sum of eigenfunctions of the form (8.7), we could establish stability if we could show that Re(m) < 0 for each of the series' terms and that the integral fo xk(x)dx exists. However, in general we cannot do so and are limited to proving stability for a restricted set of perturbations, in special cases. Note that by taking x = 0 in Equation (8.8), we obtain
m~a--~][k(O) ]=_2k(0)"
(8.13)
For cases where l / a ( 0 ) ~ 0, this automatically indicates stability for all perturbations that satisfy k ( 0 ) ~ 0. We can, however, say a great deal more.
STABILITY FOR a(x) = 1 When a(x) - 1, Equation (8.12) becomes
8.2.
(x+m+2)k"+(Ax+Am+4)k'+3Ak=O.
(8.14)
For simplicity, we rescale the problem in such a way that A = 1 and substitute z = x + m + 2 to obtain
zk" + ( z + 2)k' + 3k=O.
(8.15)
260
SHAY G U E R O N A N D SIMON A. LEVIN
Equation (8.15) can be solved analytically [by substituting k =
e-Zu(z), which transforms it to a Laguerre-type ODE]. Its two independent solutions are
k l ( z ) = e-Z(z - 2 )
(8.16)
1 1 k2(z ) = - ~ + ~-~ + e-Z(½zEi(z) - E i ( z ) ) ,
(8.17)
and
where Ei(z) is the exponential integral function. It is easy to verify that k2(z) behaves asymptotically as z -3 at z --*~. Consequently, the integral foxk(x) [see also Equation (8.3)] exists for all solutions of Equation (8.15) and for any initial conditions. This allows us to integrate both sides of Equation (8.14) from 0 to and to repeat the procedure after multiplying both sides by x. We denote the zeroth-order moment of a solution of Equation (8.14) by F 0 and its first-order moment by F 1. The initial conditions for Equation (8.14) at x = 0 are denoted by k(x = O) = k o and k ' ( x = 0) - k 0. The following identities hold: _
F o = ( m + 2 ) k ' o + ( m + 3 ) k o,
F,=(m+2)(Fo-ko).
t
(8.18)
From Equation (8.13) it follows that if m ~ - 2 we must assume that k 0 = 0. From relations (8.18) we then conclude that the conservation condition F 1 = 0 cannot be satisfied unless m = - 2 . Therefore, the only possible eigenvalue is m = - 2. It is easy to verify that, indeed, the associated eigenfunction, k(x) = e - ( x + m + 2)(x + m), satisfies the conservation condition (8.3). In this case our approach provides only limited information: Stability is ensured for initial perturbations that belong to the one-parameter family spanned by k ( x ) = e (x+m+2)(X "}'-m), but not generally. 8.3.
STABILITY
FOR a(x)= 1/x
For the case a ( x ) - 1/x, our analysis provides a better solution, although not a complete one. Equation (8.12) becomes 2 ((m*l)x+~-)k"+
[A(m+l)x+2m*4]k'+(m+3)Ak=0. (8.19)
We rescale the problem (i.e., the total population) in such a way that A = 1. Assuming m 4: - 1, we divide by m + 1, denote A = 2 / ( m + 1),
DYNAMICS OF GROUP FORMATION
261
and write Equation (8.19) as (x+A)k"+(x+A+2)k'+(l+A)k=O.
(8.20)
The asymptotic behavior of the solutions of Equation (8.20) at x --* is readily investigated by means of the W.K.B. method. One independent solution decays to 0 algebraically (asymptotically x -(A+I)-x (m+3)/(m+ 1), and the asymptotic decay of the other solution is exponential (x A- t e x ) . From Equation (8.13) we conclude that the solution of Equation (8.20) has to satisfy k(0) = 0. To prove stability for this case, we first rule out all values of m in the right half-plane as possible eigenvalues. For m such that R e ( m ) > 0, it is easy to see that 0 < R e ( A ) < 2 and that IAI<2. With no loss of generality, we assume that I m ( A ) > 0. In this open quarter circle, any solution of Equation (8.20) decays to 0 at x ~ oo. In this region, at any desired grid resolution (we used spacing of 0.01), we can investigate the solutions of Equation (8.20) numerically, with the initial conditions k(0) = O,
k'(0) = 1.
(8.21)
To avoid numerical inaccuracies that might result from the existence of dominant and subdominant "branches," we used a stiff ODE solver. Our numerical investigation reveals that the solutions behave algebraically at infinity. Since the asymptotic decay to zero is like x -(A+ 1), condition (8.3) is not satisfied (the integral does not exist). It follows that there are no eigenvalues m with nonnegative real part. For real values of m > 0, we prove the above statement analytically, as outlined below. We denote the zeroth-, first-, and second-order moments of a solution of Equation (8.20) with initial conditions (8.21) by F 0, F 1, F2, respectively. Assuming that these exist [otherwise condition (8.3) cannot be satisfied], we obtain, by integrating both sides of Equation (8.20) (after multiplying by the appropriate power of x) and using the identities (5.1) and (5.2), F°=I'
F1
A A-I'
( A - 2 ) F 2 = 0"
(8.22)
This implies that condition (8.3) is automatically satisfied for any (complex) A 4= 2 for which the moments F0, F1, F2 exist. Now, for the analytic proof that real values of m > 0 are not eigenvalues, we note that m ~ [0,m) ~ A ~ [0,21
(8.23)
262
SHAY G U E R O N AND SIMON A. LEVIN
and write Equation (8.20) in the self-adjoint form (suggested by N. Liron) 4
x+A
k=O.
(8.24)
k=0,
(8.25)
For A = 2, Equation (8.24) becomes 4
2)
x+2
and its relevant solution is k( x) = xe-X.
(8.26)
For A ~ [0, 2], we can "compare" the self-adjoint problems (8.24) and (8.25) with initial conditions (8.21) by means of Picone's theorem. It follows that between any two successive zeros of the solution of problem (8.24) there exists a zero of the solution of (8.25). Since the solution (8.26) is positive for x > 0, we conclude that the solutions of problem (8.24) are also positive for x > 0, for all A ~ [0,2]. This contradicts Equation (8.22) and therefore the assumption that F 2 exists. It follows that for A ~ [0,2] (i.e., m > 0) the solutions of the ODE (8.20) with initial conditions (8.21) behave algebraically at x ~ ~. Unfortunately, since we cannot show that the eigenfunctions are complete, we settle for showing that our problem has a continuum of eigenvalues. Of course, we do not assert that we find all of the possible eigenvalues, which is not necessary for our argument. Since F 2 exists for m ~ ( - 1,0) (A > 2), one can easily see that these are all legal eigenvalues.
To conclude, we also show that, in some cases, eigenfunctions must decay to zero monotonically. We consider the case m = - 1 that was left out above, for which Equation (8.19) takes the form k " + k ' + ~ = o.
(8.27)
The solution of Equation (8.27) that satisfies the initial conditions (8.21) is
k( x) =sin(V~x)
(8.28)
Although this solution satisfies condition (8.3), it is not a legal eigenfunction since, due to the oscillations, it violates condition (8.4).
DYNAMICS OF GROUP FORMATION 9.
263
SUMMARY
The problem of what determines group size distribution in nature is a fundamental one, carrying with it both ecological and evolutionary significance. As individual behaviors change in response to cues, so too will the distribution of groups formed. Ultimately, understanding what controls distributions is one of sorting out responses of individuals to environmental cues as well as responses to each other. Whatever the causes, as long as group sizes can be clearly defined, the fundamental processes may be captured in rates of fusion (amalgamation) and fission (fragmentation). In some cases, groups are real biological units, clearly defined and having evolutionary importance. In others, what is called a group is simply an inhomogeneity in a continuous distribution, defined in terms of a perceptual filter. Whatever the reality, the formalism applies. Paralleling methods that have proved useful in the study of polymer dynamics (e.g., [1, 4, 5, 13, 18, 21]), we develop a partial differential integral equation structured on two fundamental kernels for fusion and fission. Considering various special cases, we establish in some cases the existence and uniqueness of a steady-state distribution. Under limited conditions, the distribution may be given explicitly. The case of stability remains unresolved, although we show progress in that direction. We are pleased to acknowledge insightful discussions with Dan Cohen, Odo Diekmann, Nadav Liron, and Dan Rubenstein and the support of the Office of Naval Research through its University Research Initiative Program Award to Woods Hole Oceanographic Institute. This paper is dedicated to our late friend and colleague, Stavros Busenberg. REFERENCES 1 M. Aizenman and T. A. Bak, Convergence to equilibrium in a system of reacting polymers, Commun. Math. Phys. 65:203-230 (1979). 2 J. J. Anderson, A stochastic model for the size of fish schools, Fishery Bull. 79:315-232 (1981). 3 J.E. Cohen, Casual Groups of Monkeys and Men--Stochastic Models of Elemental Social Systems, Harvard Univ. Press, Cambridge, MA, 1971. 4 R.L. Drake, A general mathematical survey of the coagulation equation, Top. Curr. Aerosol Res. Part 2 3:201-376 (1972). 5 M.H. Ernst, Kinetics of clustering in irreversible aggregation, in Fractals in Physics: Proceedingof the Sixth International Symposium on Fractals in Physics, ICTP. Trieste, Italy, July 9-12, 1985, L. Pietronero and E. Tossati, Eds., North-Holland, Amsterdam, New York, 1986, pp. 289-302. 6 D.M. Gordon, Group-level exploration tactics in fire ants, Behavior 104:162-175 (1984).
264
SHAY GUERON AND SIMON A. LEVIN
7 D. Griinbaum, Translating stochastic density-dependent individual behavior with sensory constraints to a Eulerian model of animal swarming, J. Math. Biol. 33:139-161 (1994). 8 D. Griinbaum and A. Okubo, Modelling social animal aggregations, in Frontiers in Mathematical Ecology, Lect. Notes in Biomath Vol. 100, S. A. Levin, Ed., Springer-Verlag, Heidelberg, 1994, pp. 296-326. 9 S. Gueron and N. Liron, A model of herd grazing as a travelling wave, chemotaxis and stability, J. Math. Biol. 27:595-608 (1989). 10 E. M. Hendriks, M. H. Ernst, and R. M. Ziff, Coagulation equations with gelation, J. Stat. Phys. 31:519-563 (1983). 11 S.A. Levin, Grouping in population models, in Epidemic Models: Their Structure and Relation to Data, D. Mollison et al., Eds., Cambridge Univ. Press, 1995 (to appear). 12 S.A. Levin, A. Morin, and T. H. Powell, Patterns and process in the distribution and abundance of Antarctic krill, in Scientific Committee for the Conservation of Antarctic Marine Living Resources, Selected Scientific Papers, Part I, SCCAMLR-SSP/5, CCAMLR, Hobart, Tasmania, Australia, 1989, pp. 281-299. 13 J.H. McLeod, On an infinite set of non-linear differential equations I, Quart. J. Math. Oxford 13:119-128 (1962). 14 K.S. Norris and C. R. Schilt, Cooperative societies in three-dimensional space. On the origin of aggregations, flocks and schools with special reference to dolphins and fish, Ethol. Sociobiol. 9:149-179 (1987). 15 A. Okubo, Dynamical aspects of animal grouping: swarms, schoots, flocks and herds, Adv. Biophys. 22:1-94 (1986). 16 D. I. Rubenstein, On predation, competition, and the advantages of group living, Perspect. Ethol. 3:205-232 (1978). 17 A. R. E. Sinclair, The African Buffalo--A Study of Resource Limitation of Population, Univ. Chicago Press, Chicago, IL, 1977. 18 J.L. Spouge, An existence theorem for the discrete coagulation-fragmentation equations, Math. Proc. Camb. Phil. Soc. 96:351-356 (1984). 19 P . G . J . van Dongen and M. H. Ernst, Kinetics of reversible polymerization, J. Stat. Phys. 37:301-324 (1984). 20 P . G . V . van Dongen and M. H. Ernst, Cluster size distribution in irreversible aggregation at large times, J. Phys. A: Math. Gen. 18:2779-2793 (1985). 21 W. H. White, A global existence theorem for Smoluchowski's coagulation equations, Proc. Am. Math. Soc. 80:273-276 (1980). 22 N. G. Wolf, Behavioral ecology of herbivorous reef fishes in mixed-species foraging groups, Ph.D. Thesis, Cornell Univ., Ithaca, NY, 1983.
The Dynamics of Group Formation SHAY GUERON Department of Mathematics, Technion-- lsrael Institute of Technology, Haifa, 32000 Israel AND
SIMON A. LEVIN Department of Ecology and Evolutionary Biology, Princeton University, Princeton, New Jersey 08544-1003 Received 6 January 1994; revised 10 April 1994
ABSTRACT A general continuous model is presented for animal group size distribution. Attention is restricted to a fixed size population divided into groups of various dynamic sizes, but the approach extends easily to populations of variable size. The basic idea is to relate group size distribution to two functions, the (densitydependent) rates of fusion and fission. These functions can be estimated from data and can ultimately be related to the behavior of individuals and the dynamics of groups. For various functional forms, the stationary distributions of group sizes are sought. In several prototype cases, the stationary distribution has a peak value, the "most frequent group size," which emerges endogenously from the dynamics. The authors determine when such a peak emerges and more generally show the existence and uniqueness of the stationary distribution. Stability of stationary solutions is discussed. Progress is shown, but a general treatment remains refractory.
1.
INTRODUCTION
The problem of animal grouping is one of the most fundamental ones in biology. The nonuniform aggregation patterns of organisms ranging from small invertebrates (see, e.g., [6, 12]) to large vertebrates (e.g., [14, 16, 17, 22]) have both ecological and evolutionary significance, and the tendency to aggregate is under strong evolutionary control. The causes and consequences of grouping represent fascinating subjects for theory and will be discussed elsewhere (see, e.g., [8, 15]). In this paper, we restrict attention to relating observed grouping patterns to the processes of aggregation (fusion) and splitting (fission). In theory, both of these processes can be estimated from data and ultimately can be related to individual behavior. For example (see [11]), the decision of MATHEMATICAL BIOSCIENCES 128:243-264 (1995) © Elsevier Science Inc., 1995 0025-5564/95/$9.50 655 Avenue of the Americas, New York, NY 10010 SSDI 0025-5564(94)00074-A
244
SHAY GUERON AND SIMON A. LEVIN
an individual to join or not join a group of size n, as well as the decision by individuals in a group of size n to permit or not permit an additional individual to join can be represented in the fusion rate of groups of size 1 and groups of size n. Even the fundamental problem of pair formation/disassociation is a special case of a fission/fusion process, albeit one for which a two-sex model would be more appropriate. Some models for animal aggregation suggest that the underlying grouping mechanism involves individuals "seeking target density" (STD). In such models the individual is assumed to prefer a target number of neighbors within its perception/detection range and responds to its perceived environment accordingly. For example, based on STD, GriJnbaum [7] derives an integrodifferential equation for swarming. In a complementary approach, Gueron and Liron [9] discuss a situation where movement is influenced by a combination of chemotactic attraction and social behavior. They show that STD does not yield grouping, at least not in the restrictive sense of traveling wave solutions, whereas other social behaviors can. STD and other models that rely on local cues do not deal with group size distributions explicitly. A central question is how group size distributions in nature arise from the rules individuals follow and the rates for fusion and fission. Most important, the modal group size may be an epiphenomenon, not one preferred or even perceived by individuals. The approach presented here is not entirely novel, either within the framework of population biology or elsewhere. Cohen [3] has developed a fundamental treatment of the group size distribution problem but restricts transformations of groups to those that cause the group to increase or decrease by a single individual. Other approaches include that of Anderson [2], who presents and analyzes a linear model for the dynamics of change. Finally, a sophisticated literature exists within the study of polymers (e.g., [1, 4, 5, 10, 18-21]), in which models very similar to those we analyze herein are discussed. 2.
DEVELOPMENT OF THE BASIC MODEL
We regard the continuous model presented here as a limit of a discrete model (see, e.g., [11]). The continuous treatment facilitates calculation and introduces a number of delicate points of interpretation. For example, group sizes are allowed to vary over [0,~) (but the distribution is constrained to have an exponentially decaying "tail") although the population is of a fixed finite size. In the remainder of this section we eschew delicacy and restrict attention to the continuous model. We consider a population of total size P, distributed into groups of size x (x~[0,~)), where f(x,t) denotes the density distribution of
DYNAMICS OF G R O U P F O R M A T I O N
245
groups of size x at time t. For convenience, we drop the explicit notation for the time dependence. Instead, we distinguish time dependence by a tilde [e.g., f(x) stands for f(x, t)]. Assuming that all of the moments (about 0) of f(x) exist, we denote
Mn = fo z"f( z)dz =
the nth moment about O.
We introduce a special notation for the first two moments, since these have important biological interpretation: ~c
(~ P
=
M0
= fo f( z)dz =
=/~1 =fo
the number of groups,
zf( z)dz = the total population size.
Recall that P is fixed in our model, although it would be an easy extension to relax this assumption (see [11]). In the following sections we discuss the stationary distributions [denoted f(x)] and their associated moments (Mn). These are, by definition, time-independent and hence the tilde is omitted. The model rests on two density-dependent functions, qJ and ~b, describing rates of fusion and fission of groups. In particular,
tO(x,y) =
the probability per unit time for two groups of sizes x and y to merge. ~h(x, y ) = the probability per time for a group of size x to fragment and generate groups of size y and x - y (we ignore breakage into multiple fragments). We also introduce
p(x) =
the probability per unit time for a group of size x to break into two fragments.
Note that p(x) is not independent of the other functions but is simply (because of double counting) half the integral of 4, over all possible group sizes smaller than x [see (2.4) below]. Clearly, we can also introduce a function Q to satisfy
qb(x,y) =p(x)Q(x,y),
(2.1)
where Q(x, y) is the probability that, given that a group of size x breaks into two fragments, one of them is of size y (and the other is of size
x-y).
246
SHAY GUERON AND SIMON A. LEVIN
The quantities fb(x,y), qJ(x,y), p(x), and Q(x,y) must satisfy the following requirements:
(2.2)
(by definition),
x
(2.3)
(fragments are smaller than the "mother" group),
p(x) =3
(2.4)
(by definition), so that
foQ(X, z)dz = 2
(2.5)
qJ( x,z) = qJ( z,x)
(2.6)
[provided p(x) 4: 0] and
(symmetry). 3. A GENERAL INTEGRODIFFERENTIAL BALANCE EQUATION AND ITS STATIONARY DISTRIBUTION 3.1. TIME-DEPENDENT INTEGRODIFFERENTIAL BALANCE EQUATION The time evolution of the density distribution of groups is the result of the balance between fission (fragmentation of groups into smaller groups) and fusion (merger of groups into larger groups). It can be described by the following balance equation (see also [1]):
Ot
p(x)f(x)- fof(X)f(z)
1£.
+ -~
f(y)f(x-y)O(y,x-y)dy+ III
jxf ( y ) 4 ) ( y , x ) d y . IV
(3.1)
The four terms on the right-hand side of Equation (3.1) account for (I) fragmentation of groups of size x, (II) amalgamation of groups of size x
247
DYNAMICS OF GROUP FORMATION
with other groups, (III) generation of groups of size x by merging of two smaller groups; and (IV) fragmentation of larger groups into fragments, one of which is of size x. 3.2. THE EQUATION FOR THE STATIONARY DISTR1BUTION
A stationary distribution is a solution of Equation (3.1) that satisfies
Of(x) 3t
O.
(3.2)
The equation for this stationary distribution, f(x) (tilde omitted), is
- p(x)f(x) +-il
- f ( x ) foxy(z)
z) dz
f; f ( y ) f ( x - y ) O ( y , x - y ) d y +
ix~f ( y ) 4 , ( y , x ) d y = O .
(3.3)
We investigate the existence, uniqueness, and stability of such a stationary distribution. Furthermore, to allow comparison with data, we solve for this distribution in special cases. 4.
A MODEL TEMPLATE AND T H R E E TYPICAL EXAMPLES
For the remainder of this paper, we assume that Q(x, y) is independent of y; that is, the conditional distribution for fragments is uniform: 2
Q(x,y) =x'
(4.1)
which implies 2
ch( x,y) = x P ( X ) .
(4.2)
This somewhat restrictive assumption represents the most reasonable null case, but it will be important to consider other cases in future work. Our model is restricted to cases where the rates at which groups undergo fusion and fission depend only on the groups' sizes. Therefore, the (time-independent) functions O(x,y) and p(x) reflect the "social dynamics" of the population. We further assume that the symmetric function O(x, y) is separable, so that
qJ( x, y ) = aa( x )a( y ), where a(x) is some nonnegative function.
(4.3)
248
SHAY GUERON AND SIMON A. LEVIN
The resulting equations can be investigated numerically, even though they are resistant to analytical investigation. It is instructive, however, to consider the special case
p(x) = ~Sxa(x),
(4.4)
for which an explicit solution is feasible. As shown in the next section, this case carries substantial biological interest. Together with Equation (4.2), Equation (4.4) implies
~b(x,y) = 2/3a(x).
(4.5)
With this, the time evolution equation (3.1) becomes 0f(x,t)
Ot
Bxa( x)f( x ) - aa( x)f( x) rosa( z)f( z)dz + ½ a ~ a ( y ) f ( y ) a ( x - y ) f ( x - y)dy +2#~a(y)f(y)dy,
(4.6)
whereas the equation for the stationary distribution (3.3) translates to
- ~ x a ( x ) f ( x ) - aa(x)f(X)fo~a(z)f(z)dz + ½afo a ( y ) f ( y ) a ( x - y ) f ( x - y)dy + 2~
Jo"a ( z ) f ( z ) d z = O.
(4.7)
At this point, any specific case of the model is defined by the choice of the function a(x) and the parameters a and /3. 4.1. THREE 7)'P1CAL EXAMPLES We focus our attention on three typical examples, where the function
a(x) is decreasing, increasing, or constant. As shown below, these yield three different types of stationary distributions. Of course, other variants may easily be adopted.
Case 1 a(x) = 1.
(4.8)
In this case group fusion is independent of group size; that is, the probability for a group of size x to merge with a group of size y is
D Y N A M I C S OF G R O U P F O R M A T I O N
249
constant (independent of x and y). Furthermore, group fission [represented by p(x)] is proportional to group size, again quite a reasonable null hypothesis. Consequently, ~0(x,y)= ~, p ( x ) = / 3 x , and ~b(x,y)= 2/3 (c~ and /3 are positive constants). Case 2
(4.9)
a ( x ) = x,
which implies that ~(x, y) = axy, (o(x, y ) = 2/3x, and p ( x ) = fix 2, This accounts for grouping behavior in which large groups are more "attractive" or likely to collide and merge with others than are smaller ones, in proportion to their size. Furthermore, larger and larger groups become increasingly unstable, so that the fragmentation rate increases nonlinearly with group size. Case 3
(4.10)
a(x) = 1/x,
which implies that ~b(x,y)= a / x y , ¢ ( x , y ) = 2 / 3 / x , and p ( x ) = / 3 . Thus, groups break at a rate independent of their size, for example, due to external disturbances that affect all groups uniformly. Furthermore, the functional form of ~O means that solitary individuals, or individuals in small groups, do not remain in such groups for long. Possibly, due to the fact that it is disadvantageous to be in small groups, small groups merge at high rates to form larger ones. 5.
IDENTITIES, R E C U R R E N C E RELATIONS, AND CONVERGENCE PROPERTIES
5.1. IDENTITIES Here we prove two identities that are used later. Let f ( x ) be a function, defined for x~[0,~), that decays to 0 exponentially as x--* ~. Then all of its moments about zero (denoted M,) exist, and we obtain
=L
I(y)dyf,
+ yl"f(z)dz n
=
C.MkM,,_ k k=0
(5.1/
250
SHAY GUERON
A N D S I M O N A. L E V I N
and
fo xn d~ fx f ( z ) dz = - n +1l 5.2.
( z ) dz Sozx" d~ fo~z " + l f ( z ) d z = n +1l M,+ 1. (5.2)
DERIVING R E C U R R E N C E R E L A T I O N S FOR THE M O M E N T S
We define the auxiliary function
(5.3)
e(*) = a ( x ) h x ) and denote its moments about zero by Nn = f0~xn~(x) dx
(5.4)
(assuming that N,, M,, N,, and M n exist for all natural n and for all t). A recurrence relation between the moments of f ( x ) and g(x) can be obtained by multiplying both sides of Equation (4.6) by x ~ and integrating from 0 to ~ with respect to x. Applying the identities (5.1) and (5.2), we obtain
d& .
dt
.
- 1., + lf. l ~n N
aNolqn +½a ~
C nkN" k N "n _ k
k=O n-1
=-
n - 1 N# - + i +½a E /3 ~---~
k" "
CnNkNn_ k .
(5.5)
k=l
For n = 0, one must use only the "explicit" (first) equality in Equation (5.5). For n = 1, Equation (5.5) is an identity with both sides equal to zero. This is a manifestation of the conservation of the total population size (/~a). Note that for our three special cases defined above, we have /Vn = h4n
[for a(x) = 1],
(5.6a)
N~ =/~n÷l
[for a(x) = x],
(5.6b)
[for a(x) = 1/x].
(5.6c)
N~÷1 = / ~ 5.3.
C O N V E R G E N C E P R O P E R T Y FOR THE CASE a ( x ) = 1
The case a(x)= 1 is the only case where the system of coupled differential equations (5.5) yields a straightforwardly solvable equation, at least with our skills. We derive an equation for the evolution (in
DYNAMICS OF GROUP FORMATION
251
time) of the number of groups, (~ ( = M0), by using the identity (5.5) for n -- 0. From (5.6) we have M n -- ]Qn, which gives dG
1
dt = - T A G
- 2
+tiP
(P=2Q,).
(5.7)
The solution of Equation (5.7) can be found analytically by substitution of
U
=½.d.
(5.8)
This yields
mexp{[(1/2)tap]'/zt}_nexp{_[(1/2)~ap]'/2t} (~(t) = i -2t-Pa m e x p { [ ( 1 / 2 ) f l a P ] '/2 t} + n exp{- [ ( 1 / 2 ) f l a P ] 1/2t} (5.9) where the ratio m / n is determined by the initial conditions. It is now easy to see, as was originally shown by Aizenman and Bak [1], that regardless of initial conditions the number of groups converges to the limit lim (~(t) = i 2fl--e . 6.
(5.10)
EXISTENCE AND UNIQUENESS OF THE STATIONARY DISTRIBUTION
Although we cannot solve the more general problem or prove stability, we can establish the existence and uniqueness of the stationary distribution for Equations (4.1)-(4.7). Note that these are more general than the three special cases mentioned earlier. To that end, it is important to notice first that if one defines the auxiliary function g(x) by
g(x) - a ( x ) f ( x )
(6.1)
and writes Equation (4.7) in terms of g(x), it is reduced to the form of Equation (4.7) written in terms of f(x) for the case a(x) -- 1. Therefore, we conclude that showing the existence and uniqueness of the stationary distribution for the special case a(x) = 1 implies the same result for the more general case.
252
SHAY GUERON AND SIMON A. LEVIN
We denote the Laplace transform of f(x) by
fo e-SZf( z)dz - F( s).
(6.2)
F ( 0 ) = G.
(6.3)
zc
Sa( f( x) ) = For s = 0 we have
To solve the equation for the stationary distribution we take the Laplace transform of Equation (4.7) and obtain
/3F'(s)- o~GF(s)+½o~FZ(s)+2/3G - 2 / 3 F ( s ) = 0 . s
(6.4)
s
Equation (6.4) is a first-order O D E with a singular point at s = 0. For every s o > 0, this equation is regular on the interval [s0,~). As such, all of its solutions on this interval belong to a one-parameter family, namely, /3 1 F ( s ) = 2 S (s-~--d).
(6.5)
Continuing a solution from the above family to s ~ [0, ~) determines the parameter d uniquely, and it is easy to verify that the only solution of Equation (6.4) on the interval s ~ [0, ~) is s +2(/3/o~G)
"
(6.6)
Transforming F(s) back to the real domain, we obtain the stationary distribution, which is uniquely determined by G:
f( x) = 2 ~ e -2(tJ/'~c)x OL
(6.7)
Now, instead of fixing G, one can equivalently choose P as the normalization factor. For a given value of P, the stationary solution is
and the number of groups (G) is related to P by
G = ~/-~~
(6.9)
DYNAMICS OF GROUP FORMATION
253
Returning to the general case we write the stationary solution e
where the value of h is determined by the equation P=2--~
a - z~
e AZdz.
(6.11)
This result is valid for any positive a(x) with algebraic asymptotic behavior at x ~ ~ (to ensure that all moments exist). Singularity occurs when a(0)= 0 [e.g., for a(x)= x], although the integral may still converge. Consequently, to ensure that the population (P) remains finite, we must restrict the choice of a(x) to satisfy lim0 -a(7 -x ) > 0. ,-~
(6.12)
All three special cases satisfy these conditions. Figure 1 displays the shapes of the relationships (6.10) for the three special cases, with a, /3, and P all set equal to 1. 7.
I N T E R P R E T A T I O N OF T H E STATIONARY SOLUTIONS
From Equation (6.10) it is easy to see that the properties of the resulting stationary solution depend only on a(x). Also, it is easy to integrate the expressions for f(x) and xf(x) over any desired interval [xl,x 2] to obtain the number of groups whose sizes are in this interval and the associated total population. Equation (6.10) demonstrates that for every choice of a nondecreasing function a(x) [e.g., a(x)= 1 and a(x)= x] the resulting stationary distribution is a monotonically decreasing function of x. For a(x)= 1 we obtain an exponential distribution of group sizes in which the most frequent groups are the "small ones." The same conclusion is true for a(x)= x, where the resulting distribution f(x) is singular at x = 0. In this case we have small groups piling up near x = 0 and larger groups becoming less abundant. On the other hand, we conclude that any choice of a decreasing function a(x) [e.g., a(x) = 1/x] yields a group size stationary distribution that has at least one local maximum. This maximum is the modal group size. It follows, not surprisingly, that a population can display a modal group size without having this quantity built explicitly into the dynamics.
254
SHAY GUERON AND SIMON A. LEVIN
2.5-
0
0.5
1
1.5
2
2.5
3
3.5
(a~ FIG. 1. The shapes of the stationary solutions in (6.10) for the three special cases (a) a ( x ) = 1, (b) a ( x ) = x, (c) a(x)= 1/x. In all three plots, a, /3, and P are set equal to 1. The resulting values of h are (a) h. = Vr2, (b) h. = 2, (c) h = 41/3.
DYNAMICS OF GROUP FORMATION
ol 0
I 0.5
255
t
t
1
1.5
x
(b) FIG. 1. Continued.
256
SHAY GUERON
A N D S I M O N A. L E V I N
iT
0.8-
0.6-
0.4-
,
\ 0.2-
\
\
\\ \
\
h
2 x
(c) FIG. 1. Continued.
L
3
t
4
DYNAMICS OF GROUP FORMATION 8.
257
STABILITY ANALYSIS
In this section we present work in progress; the full problem remains impenetrable. We investigate the local stability of the stationary solutions for the cases a(x) = 1 and a(x) = x. As demonstrated below, this leads to a nontrivial eigenvalue problem, for which we have only a partial solution. Our hope is that others will be inspired to complete the proof, and hence we present these results as open questions. 8.1. EIGENVALUE PROBLEM FOR A GENERAL PERTURBATION We start with the dynamical balance PDE (3.1) written in terms of ff,(x) - f(x)a(x). With no loss of generality we assume ~ =/3 = 1 (this can be achieved by time and group size rescaling) and obtain
~,(x) a( x)
x~,( x)
'S/
+g We consider 2e-AX/a(x),
-
~,(x)f° g(y)dy
g(x-y)g(y)dy+2
a perturbation
f(x) =
on the
/;
g(y)dy.
stationary
2e :ax + h ( x ) a(x) '
solution
(8.1) f(x)=
(8.2)
for some (time-dependent) function h(x). Conservation of the population size restricts the initial perturbation /~(x)t=0, and the integral constraint
t
X~a(x)],=o
=o
(8.3)
must hold. To allow for linearization, we start with /Tt(x),=0 <
(8.4)
258
SHAY GUERON AND SIMON A. LEVIN
We linearize Equation (8.3) by omitting the two quadratic terms
fo~(x-y)~(y)dy,
[~(X)fo~h(y)dy
(8.5)
and obtain the linear equation, which, after reorganization, can be written as
-
+ 2e-AX~t(y)e ~' dy + 2fx~[t(y)dy.
(8.6)
Assuming that/~ can be expanded in terms of eigenfunctions, we use the standard separation of variables
tl = emtk(x)
(8.7)
and obtain an equation for k(x):
2 m k(x) a(x) =_ xk(x)--xk
_2e_.Xfok(y)dy
+ 2e-~Xfok(y)eXY dy + 2fx k(y)dy.
(8.8)
Defining the linear operator 5r by k ( x ) ) - - x k ( x ) - -2 2k
_2e_~fo~k(y)dy ~z
+2e-~Xfok(y)eXYdy+2fx k(y)dy,
(8.9)
we regard the problem of solving Equation (8.8) as being one of finding "weighted" eigenvalues of the operator 5r,
~( k( x) ) = m k( x) a(x) [where the weight function is complex.
1/a(x)].
(8.10)
Note that m is in general
DYNAMICS OF GROUP FORMATION
259
We point out that 5r is a "conserving operator." Multiplying both sides of Equation (8.8) by x and using the identities (5.1) and (5.2), we have ~c
fo x~(k(x))dx = 0.
(8.11)
As a result, for any eigenfunction k(x) that satisfies (8.10), the perturbation component of the form (8.7) conserves the population size automatically. The integral equation (8.8) can be converted to a second-order ODE. This is achieved by differentiating it twice with respect to x (assuming the necessary smoothness), which yields
m[a--~
+
+ Am k(x) a(x) ] + 3 A k ( x )
(Ax+4)k'(x) 0,
(8.12)
where a prime denotes differentiation with respect to x. For a given choice for the function a(x), by expanding the solution as a sum of eigenfunctions of the form (8.7), we could establish stability if we could show that Re(m) < 0 for each of the series' terms and that the integral fo xk(x)dx exists. However, in general we cannot do so and are limited to proving stability for a restricted set of perturbations, in special cases. Note that by taking x = 0 in Equation (8.8), we obtain
m~a--~][k(O) ]=_2k(0)"
(8.13)
For cases where l / a ( 0 ) ~ 0, this automatically indicates stability for all perturbations that satisfy k ( 0 ) ~ 0. We can, however, say a great deal more.
STABILITY FOR a(x) = 1 When a(x) - 1, Equation (8.12) becomes
8.2.
(x+m+2)k"+(Ax+Am+4)k'+3Ak=O.
(8.14)
For simplicity, we rescale the problem in such a way that A = 1 and substitute z = x + m + 2 to obtain
zk" + ( z + 2)k' + 3k=O.
(8.15)
260
SHAY G U E R O N A N D SIMON A. LEVIN
Equation (8.15) can be solved analytically [by substituting k =
e-Zu(z), which transforms it to a Laguerre-type ODE]. Its two independent solutions are
k l ( z ) = e-Z(z - 2 )
(8.16)
1 1 k2(z ) = - ~ + ~-~ + e-Z(½zEi(z) - E i ( z ) ) ,
(8.17)
and
where Ei(z) is the exponential integral function. It is easy to verify that k2(z) behaves asymptotically as z -3 at z --*~. Consequently, the integral foxk(x) [see also Equation (8.3)] exists for all solutions of Equation (8.15) and for any initial conditions. This allows us to integrate both sides of Equation (8.14) from 0 to and to repeat the procedure after multiplying both sides by x. We denote the zeroth-order moment of a solution of Equation (8.14) by F 0 and its first-order moment by F 1. The initial conditions for Equation (8.14) at x = 0 are denoted by k(x = O) = k o and k ' ( x = 0) - k 0. The following identities hold: _
F o = ( m + 2 ) k ' o + ( m + 3 ) k o,
F,=(m+2)(Fo-ko).
t
(8.18)
From Equation (8.13) it follows that if m ~ - 2 we must assume that k 0 = 0. From relations (8.18) we then conclude that the conservation condition F 1 = 0 cannot be satisfied unless m = - 2 . Therefore, the only possible eigenvalue is m = - 2. It is easy to verify that, indeed, the associated eigenfunction, k(x) = e - ( x + m + 2)(x + m), satisfies the conservation condition (8.3). In this case our approach provides only limited information: Stability is ensured for initial perturbations that belong to the one-parameter family spanned by k ( x ) = e (x+m+2)(X "}'-m), but not generally. 8.3.
STABILITY
FOR a(x)= 1/x
For the case a ( x ) - 1/x, our analysis provides a better solution, although not a complete one. Equation (8.12) becomes 2 ((m*l)x+~-)k"+
[A(m+l)x+2m*4]k'+(m+3)Ak=0. (8.19)
We rescale the problem (i.e., the total population) in such a way that A = 1. Assuming m 4: - 1, we divide by m + 1, denote A = 2 / ( m + 1),
DYNAMICS OF GROUP FORMATION
261
and write Equation (8.19) as (x+A)k"+(x+A+2)k'+(l+A)k=O.
(8.20)
The asymptotic behavior of the solutions of Equation (8.20) at x --* is readily investigated by means of the W.K.B. method. One independent solution decays to 0 algebraically (asymptotically x -(A+I)-x (m+3)/(m+ 1), and the asymptotic decay of the other solution is exponential (x A- t e x ) . From Equation (8.13) we conclude that the solution of Equation (8.20) has to satisfy k(0) = 0. To prove stability for this case, we first rule out all values of m in the right half-plane as possible eigenvalues. For m such that R e ( m ) > 0, it is easy to see that 0 < R e ( A ) < 2 and that IAI<2. With no loss of generality, we assume that I m ( A ) > 0. In this open quarter circle, any solution of Equation (8.20) decays to 0 at x ~ oo. In this region, at any desired grid resolution (we used spacing of 0.01), we can investigate the solutions of Equation (8.20) numerically, with the initial conditions k(0) = O,
k'(0) = 1.
(8.21)
To avoid numerical inaccuracies that might result from the existence of dominant and subdominant "branches," we used a stiff ODE solver. Our numerical investigation reveals that the solutions behave algebraically at infinity. Since the asymptotic decay to zero is like x -(A+ 1), condition (8.3) is not satisfied (the integral does not exist). It follows that there are no eigenvalues m with nonnegative real part. For real values of m > 0, we prove the above statement analytically, as outlined below. We denote the zeroth-, first-, and second-order moments of a solution of Equation (8.20) with initial conditions (8.21) by F 0, F 1, F2, respectively. Assuming that these exist [otherwise condition (8.3) cannot be satisfied], we obtain, by integrating both sides of Equation (8.20) (after multiplying by the appropriate power of x) and using the identities (5.1) and (5.2), F°=I'
F1
A A-I'
( A - 2 ) F 2 = 0"
(8.22)
This implies that condition (8.3) is automatically satisfied for any (complex) A 4= 2 for which the moments F0, F1, F2 exist. Now, for the analytic proof that real values of m > 0 are not eigenvalues, we note that m ~ [0,m) ~ A ~ [0,21
(8.23)
262
SHAY G U E R O N AND SIMON A. LEVIN
and write Equation (8.20) in the self-adjoint form (suggested by N. Liron) 4
x+A
k=O.
(8.24)
k=0,
(8.25)
For A = 2, Equation (8.24) becomes 4
2)
x+2
and its relevant solution is k( x) = xe-X.
(8.26)
For A ~ [0, 2], we can "compare" the self-adjoint problems (8.24) and (8.25) with initial conditions (8.21) by means of Picone's theorem. It follows that between any two successive zeros of the solution of problem (8.24) there exists a zero of the solution of (8.25). Since the solution (8.26) is positive for x > 0, we conclude that the solutions of problem (8.24) are also positive for x > 0, for all A ~ [0,2]. This contradicts Equation (8.22) and therefore the assumption that F 2 exists. It follows that for A ~ [0,2] (i.e., m > 0) the solutions of the ODE (8.20) with initial conditions (8.21) behave algebraically at x ~ ~. Unfortunately, since we cannot show that the eigenfunctions are complete, we settle for showing that our problem has a continuum of eigenvalues. Of course, we do not assert that we find all of the possible eigenvalues, which is not necessary for our argument. Since F 2 exists for m ~ ( - 1,0) (A > 2), one can easily see that these are all legal eigenvalues.
To conclude, we also show that, in some cases, eigenfunctions must decay to zero monotonically. We consider the case m = - 1 that was left out above, for which Equation (8.19) takes the form k " + k ' + ~ = o.
(8.27)
The solution of Equation (8.27) that satisfies the initial conditions (8.21) is
k( x) =sin(V~x)
(8.28)
Although this solution satisfies condition (8.3), it is not a legal eigenfunction since, due to the oscillations, it violates condition (8.4).
DYNAMICS OF GROUP FORMATION 9.
263
SUMMARY
The problem of what determines group size distribution in nature is a fundamental one, carrying with it both ecological and evolutionary significance. As individual behaviors change in response to cues, so too will the distribution of groups formed. Ultimately, understanding what controls distributions is one of sorting out responses of individuals to environmental cues as well as responses to each other. Whatever the causes, as long as group sizes can be clearly defined, the fundamental processes may be captured in rates of fusion (amalgamation) and fission (fragmentation). In some cases, groups are real biological units, clearly defined and having evolutionary importance. In others, what is called a group is simply an inhomogeneity in a continuous distribution, defined in terms of a perceptual filter. Whatever the reality, the formalism applies. Paralleling methods that have proved useful in the study of polymer dynamics (e.g., [1, 4, 5, 13, 18, 21]), we develop a partial differential integral equation structured on two fundamental kernels for fusion and fission. Considering various special cases, we establish in some cases the existence and uniqueness of a steady-state distribution. Under limited conditions, the distribution may be given explicitly. The case of stability remains unresolved, although we show progress in that direction. We are pleased to acknowledge insightful discussions with Dan Cohen, Odo Diekmann, Nadav Liron, and Dan Rubenstein and the support of the Office of Naval Research through its University Research Initiative Program Award to Woods Hole Oceanographic Institute. This paper is dedicated to our late friend and colleague, Stavros Busenberg. REFERENCES 1 M. Aizenman and T. A. Bak, Convergence to equilibrium in a system of reacting polymers, Commun. Math. Phys. 65:203-230 (1979). 2 J. J. Anderson, A stochastic model for the size of fish schools, Fishery Bull. 79:315-232 (1981). 3 J.E. Cohen, Casual Groups of Monkeys and Men--Stochastic Models of Elemental Social Systems, Harvard Univ. Press, Cambridge, MA, 1971. 4 R.L. Drake, A general mathematical survey of the coagulation equation, Top. Curr. Aerosol Res. Part 2 3:201-376 (1972). 5 M.H. Ernst, Kinetics of clustering in irreversible aggregation, in Fractals in Physics: Proceedingof the Sixth International Symposium on Fractals in Physics, ICTP. Trieste, Italy, July 9-12, 1985, L. Pietronero and E. Tossati, Eds., North-Holland, Amsterdam, New York, 1986, pp. 289-302. 6 D.M. Gordon, Group-level exploration tactics in fire ants, Behavior 104:162-175 (1984).
264
SHAY GUERON AND SIMON A. LEVIN
7 D. Griinbaum, Translating stochastic density-dependent individual behavior with sensory constraints to a Eulerian model of animal swarming, J. Math. Biol. 33:139-161 (1994). 8 D. Griinbaum and A. Okubo, Modelling social animal aggregations, in Frontiers in Mathematical Ecology, Lect. Notes in Biomath Vol. 100, S. A. Levin, Ed., Springer-Verlag, Heidelberg, 1994, pp. 296-326. 9 S. Gueron and N. Liron, A model of herd grazing as a travelling wave, chemotaxis and stability, J. Math. Biol. 27:595-608 (1989). 10 E. M. Hendriks, M. H. Ernst, and R. M. Ziff, Coagulation equations with gelation, J. Stat. Phys. 31:519-563 (1983). 11 S.A. Levin, Grouping in population models, in Epidemic Models: Their Structure and Relation to Data, D. Mollison et al., Eds., Cambridge Univ. Press, 1995 (to appear). 12 S.A. Levin, A. Morin, and T. H. Powell, Patterns and process in the distribution and abundance of Antarctic krill, in Scientific Committee for the Conservation of Antarctic Marine Living Resources, Selected Scientific Papers, Part I, SCCAMLR-SSP/5, CCAMLR, Hobart, Tasmania, Australia, 1989, pp. 281-299. 13 J.H. McLeod, On an infinite set of non-linear differential equations I, Quart. J. Math. Oxford 13:119-128 (1962). 14 K.S. Norris and C. R. Schilt, Cooperative societies in three-dimensional space. On the origin of aggregations, flocks and schools with special reference to dolphins and fish, Ethol. Sociobiol. 9:149-179 (1987). 15 A. Okubo, Dynamical aspects of animal grouping: swarms, schoots, flocks and herds, Adv. Biophys. 22:1-94 (1986). 16 D. I. Rubenstein, On predation, competition, and the advantages of group living, Perspect. Ethol. 3:205-232 (1978). 17 A. R. E. Sinclair, The African Buffalo--A Study of Resource Limitation of Population, Univ. Chicago Press, Chicago, IL, 1977. 18 J.L. Spouge, An existence theorem for the discrete coagulation-fragmentation equations, Math. Proc. Camb. Phil. Soc. 96:351-356 (1984). 19 P . G . J . van Dongen and M. H. Ernst, Kinetics of reversible polymerization, J. Stat. Phys. 37:301-324 (1984). 20 P . G . V . van Dongen and M. H. Ernst, Cluster size distribution in irreversible aggregation at large times, J. Phys. A: Math. Gen. 18:2779-2793 (1985). 21 W. H. White, A global existence theorem for Smoluchowski's coagulation equations, Proc. Am. Math. Soc. 80:273-276 (1980). 22 N. G. Wolf, Behavioral ecology of herbivorous reef fishes in mixed-species foraging groups, Ph.D. Thesis, Cornell Univ., Ithaca, NY, 1983.