From Ants to Fishing Vessels: A Simple Model for Herding and Exploitation of Finite Resources
FFrom Ants to Fishing Vessels: A Simple Model for Herding and Exploitation ofFinite Resources
José Moran,
1, 2, 3, 4, ∗ Antoine Fosset,
1, 3
Alan Kirman, and Michael Benzaquen
1, 3, 5 Chair of Econophysics and Complex Systems, Ecole polytechnique, 91128 Palaiseau Cedex, France Centre d’Analyse et de Mathématique Sociales, EHESS, 54 Boulevard Raspail, 75006 Paris, France LadHyX UMR CNRS 7646, Ecole polytechnique, 91128 Palaiseau Cedex, France Complexity Science Hub Vienna, Josefstädter Straße 39, A-1080, Austria Capital Fund Management, 23 Rue de l’Université, 75007 Paris, France (Dated: September 23, 2020)We analyse the dynamics of fishing vessels with different home ports in an area where these vessels, inchoosing where to fish, are influenced by their own experience in the past and by their current observationof the locations of other vessels in the fleet. Empirical data from the boats near Ancona and Pescara showsstylized statistical properties that are reminiscent of Kirman and Föllmer’s ant recruitment model, althoughwith two ant colonies represented by the two ports. From the point of view of a fisherman, the two fishingareas are not equally attractive, and he tends to prefer the one closer to where he is based. This pieceof evidence led us to extend the original ants model to a situation with two asymmetric zones and finiteresources. We show that, in the mean-field regime, our model exhibits the same properties as the empiricaldata. We obtain a phase diagram that separates high and low herding regimes, but also fish populationextinction. Our analysis may have interesting policy implications for the ecology of fishing areas.
I. INTRODUCTION
A problem of general interest is that of the individual and collective exploitation of a resource. Dependingon the particular context, the dynamics can be very different. A crucial factor is the effect of the behaviour ofindividuals on the collective outcome. In financial markets for example, the decision to buy may enhance thevalue of the resource for others as the price of an asset may increase as the demand for it grows. This positivefeedback can lead to “herd behaviour” and to creation of “bubbles”. If, on the other hand, the resource is infixed supply or can only generate a limited flow, as in the case of agricultural production, over exploitation canlead to its exhaustion when individuals do not take account of the overall consequences of their actions. Thisleads to what has been called “The Tragedy of the Commons” in [ ] .In this paper we use a version of a model which was developed in the context of financial markets but wemodify it to look at a problem of exhaustible resources, in particular that of fisheries. There is a substantialliterature on fishing management which analyses the causes of over exploitation and the behaviour that leadsto this. Much of that literature was based on understanding the strategies that individual boats use to decidewhen and where to fish. The simplest idea is that the individuals base their decisions on Catch per Unit Effort(CPUE), see [ ] . This suggests that boats fish until their catch falls below a certain threshold and then moveon. This is a purely individualistic model and argues that past individual experience is an adequate basis fordecision making. Two questions arise here. Firstly, can one deduce the aggregate behaviour from the observedbehaviour of individual vessels, and secondly, does the behaviour of other vessels influence the choices of aparticular boat? The answer to the first question lies in the development of satellite technology which allowsindividual vessels to be identified and followed; this information provides a basis for analysing the individualand collective behaviour of fishing fleets. It is, of course, known that vessels do not act in total isolation anda model using tracking data for New Zealand fisheries was, for example, studied in Ref. [ ] . This came tothe conclusion that “there is evidence that vessels make decisions about where to fish based on both theirown recent catch history and on observation about the location and aggregation of other vessels. There is no ∗ [email protected] a r X i v : . [ phy s i c s . s o c - ph ] S e p evidence that there is enough information transfer for vessels to make decisions on the basis of catch rates ofthe other vessels in the fleet”. What was suggested was that while the influence of other players is taken intoaccount, because of the limited information about the performance of other vessels it may not be the majordriving influence for collective behaviour.However, a more radical approach, abandoning a simple optimization approach had been developed earlierby Allen and McGlade [ ] . They developed models in part based on the Lotka-Volterra equations which alreadyincorporated recent advances in the understanding of the evolution of complex systems. They studied herdbehaviour and simulations of a dynamic model of a Nova Scotia fishery. Their analysis revealed that humanresponses amplify rapid random fluctuations in recruitment and excite strong Lotka-Volterra type oscillationsin a system that would normally settle to a stable stationary state. Their dynamic, multi-species, multi-fleetspatial model was calibrated to the Nova Scotian groundfish fisheries. They examined the role of “exploration”and “exploitation”. They identified two types of hunters, “stochasts” or high-risk takers, and “cartesian” fol-lowers, or low risk takers. The result of the interaction between the two reveals, as they say, “the ‘out of phase’relationship between abundance and the ease with which fishermen locate a highly sought species and its con-verse”. They emphasize, contrary to more conventional analysis, “the importance of information exchange indefining the attractivity of a particular fishing zone to different fleets and the ability of the model to take intoaccount coded information, misinformation, spying and lying; and the fact that models based on global princi-ples, such as ‘optimal efficiency’ or ‘maximum profit’, are clearly of dubious relevance to the real world.” Thecrucial difference between this and the work previously cited is that much more weight is given to informationabout the activity of others and the content of the messages about that activity is assumed to be much richer.Our approach is in this spirit and is based on a model in which agents are “recruited” to a source of profitby those already benefiting from that source. The actors follow simple rules but their interaction can produceinteresting dynamics. A related approach by computer scientists [ ] suggested that the result might be thatof a uniform distribution across the space in which the resource is found. We show that, depending on theweight given to the behaviour of others, vessels can typically operate near to their home port with occasionalexcursions to another area, but that changing the parameters of the model can lead to a persistent mixing of thetwo fleets with some boats from each area fishing in the other area. Since what is important is the probabilitythat a boat follows others, the distribution of the boats over the two areas is determined by a stochastic process.This recalls a result of Allen and McGlade in which the survival of the fishery was dependent on the existenceof some vessels which chose the place to fish at random and, as in many models of interaction, a degree ofrandomness may make an important contribution to the overall dynamics of the system.The paper is organised as follows. In Section II, we present an overview of the data used in this article andintroduce all relevant definitions. Section III introduces a model intended to reproduce the main stylized factspresent in the data. The model is developed in Section IV using simplifying assumptions that are justified bynumerical simulations. Finally, in Section V, we discuss further consequences of the model and in particularthe different collective “phases” that describe the aggregate behaviour of fishing vessels. II. EMPIRICAL FISHING DATA
As mentioned above, while applicable to a wider range of situations, our work was originally inspired byimitation and herding effects in fishing areas. Here we present the data we use together with some stylizedfacts, both quantitative and qualitative.
A. Description of the data
We use the Fishing Vessels Dataset from Global Fishing Watch [ ] from Octobre 2012 to December 2016.Since our aim is to analyse the behaviour of fishermen seeking to exploit clearly distinguishable fishing areas,we geographically focus on the Adriatic Sea and specifically on the area encompassing the Italian cities ofAncona and Pescara in which two of the largest fishing harbours and fish markets are situated (see for example [ ] for a detailed study and description of the Ancona fishing market). The two cities are separated by areasonable distance of about 150 km, meaning that boats based in one city can easily find themselves fishingclose to the other. Further, the existence of large and comparable fish markets in both cities hints at thepossibility of matching fishing activity to market data, provided of course one has access to the latter. Notethat while another city, San Benedetto del Tronto, lies between Ancona and Pescara, it is responsible for arather negligible amount of the activity in the area.We have also restricted our analysis to the behaviour of trawlers. These boats have a low cruise speed andfish in shallow waters close to the coast. A reasonable hypothesis, which we have confirmed with the localmarket authorities, is that trawlers fishing in the area are based in either one of the two cities and go out fora short amount of time before coming back to sell their catch on the local market. In particular we were toldthat, due to the policy of the market to sell fresh local fish, vessels (almost) always get back to the port after24 hours. We were also told that while there is no ban for a boat registered in a given port to land their fishelsewhere, this seldom happens. In other words, one expects trawlers based in, say, Pescara to leave port, fishfor a day and then come back to sell their catch.The reduced data set consists of daily tracking of these trawlers, identified by their 9-digit Maritime MobileService Identity (MMSI) number. Each vessel is tracked on a latitude-longitude grid with resolution 0.1.1squared degrees. At Ancona and Pescara’s latitude ( ≈ ◦ North), this implies a spatial resolution of ≈ × (latitude by longitude). Finally, a preliminary study of the data shows that there is a significant reductionof fishing activity from Friday to Sunday, consistent with markets being open Monday through Thursday only.We have thus dropped the former from our data set, keeping only trading days to ensure significant fishingactivity. B. Defining fishing areas
To assign each trawler to its base port (Ancona or Pescara), we use the following heuristic procedure, whichwe then cross-validate with MMSI data provided by the Ancona market authorities. We introduce the notations: • h i ( x , t ) the time spent by trawler i fishing at grid-point x on day t , • w i ( x ) : = (cid:80) t h i ( x , t ) / (cid:80) y , s h i ( y , s ) , for the average fraction of time spent by trawler i fishing at point x , • d A ( x ) the distance between point x and Ancona, and d the distance between the two cities, • d i A : = (cid:80) x w i ( x ) d A ( x ) , the average distance separating trawler i and Ancona when it is fishing, • D i A : = (cid:80) x w i ( x )[ d A ( x )] , the average square distance between trawler i and Ancona,and of course symmetrically for Pescara with index P. We then define the neighborhood of Ancona and Pescara asthe pseudo-ellipsoid with focal points the two ports, i.e. the set { x | d A ( x ) + d P ( x ) ≤ d } , of course excludingland, see Fig. 1(c). We restrict our analysis to trawlers evolving within this area, namely { i | D i A + D i P ≤ d } . Wethen assign the trawlers to one of the two ports according to their average distance to each of them. Definingtwo distinct areas as: (cid:68) A = (cid:8) x | d A ( x ) ≤ d P ( x ) and d A ( x ) + d P ( x ) ≤ d (cid:9) (1a) (cid:68) P = (cid:8) x | d P ( x ) < d A ( x ) and d A ( x ) + d P ( x ) ≤ d (cid:9) , (1b)a given trawler is assigned to, say, Pescara if its fishing time-weighted average position lies in (cid:68) P . In other words i ∈ Pescara (resp. Ancona) if d i P ≤ d i A (resp. d i A < d i P ) . To validate our method of home port identification, we Note that with this publicly available database, our analysis can be reproduced in any other place of the world where two competingharbours lie reasonably close to each other. According to the director of the Ancona fish market, there are no relationships with nearby wholesale markets (Pescara and SanBenedetto del Tronto), and, two or three times a year, it happens that a boat based in the nearby port in the north (Fano or Cattolica)comes to sell. (a) n A − n P (b) N A + N P (c) τ (days) − − C − n A ( τ ) (d) Ancona Exponential fitEmpirical τ (days) − − C − n P ( τ ) (e) Pescara Exponential fitEmpirical
FIG. 1. Description of empirical data. Blue curves and markers correspond to data related to the area of Ancona, while redcurves and markers correspond to Pescara. (a) Plot of the fractions n i ( t ) , as defined in Eq. (2) (b) Plot of the total numberof active boats through time N A + N P . (c) Satellite view of the Adriatic Sea along with the areas (cid:68) A and (cid:68) P , as definedin Eq. (1). (d) and (e) Autocorrelation functions C − n i ( τ ) , as defined in Eq. (3) for both zones. For Ancona we find anexponential fit with a decay rate of ≈
11 days, while for Pescara we find a decay of ≈
33 days. were able confront our classification to the list of the Ancona-based trawlers, kindly provided by the Anconafish market authorities. Up to a few minor errors, notably related to having identified as Ancona-based a fewvessels based in the much smaller San Benedetto del Tronto, the cross-check was successful. Over the wholeperiod we counted N A =
108 Ancona-based trawlers and N P =
118 Pescara-based trawlers.
C. Stylized facts
Having tagged each boat to either Ancona or Pescara, we now turn to studying the dynamics of fishing withinthe two areas (cid:68) A and (cid:68) P . We define the fraction n A ( t ) of time spent by Ancona-based vessels fishing in (cid:68) A namely: n A ( t ) = (cid:80) x ∈(cid:68) A , i ∈ Ancona h i ( x , t ) (cid:80) y , i ∈ Ancona h i ( y , t ) , (2)and vice-versa n P ( t ) for Pescara. Figure 1(b) displays the evolution of n A ( t ) and n B ( t ) throughout the periodof interest. While these fractions are most often very close to 1, indicating as one would intuitively expect n P < ( n ) AnconaBeta distribution fitEmpirical n Pescara
FIG. 2. Cumulative distribution function (cdf) of the fractions n A and n P as defined in Eq. (2). The solid red curvescorrespond to a fit with a generalized Beta distribution, which has a cdf given by P > ( n ) = C (cid:82) n d x x γ − ( − x ) γ − with C anormalization constant. The parameters for Ancona are γ = γ = γ = γ = that trawlers spend most of their time fishing near their home port, one can see, however, that they regularlyundergo persistent excursions, revealing that a sizeable fraction of the vessels in each area decide collectivelyto go elsewhere.To evaluate the typical length of such excursions, Figs. 1 (d) and (e) display the auto-correlation functions: C − n ( τ ) : = Cor ( − n ( t + τ ) , 1 − n ( t )) , (3)for both n A ( t ) and n P ( t ) . These are well fitted by the sum of a delta-peak at 0, which can be attributed tomeasurement noise and other exogenous factors such as the weather, and an exponentially decaying functionwith typical timescale ranging from ≈
11 to ≈
30 days. Interestingly enough, Fig. 2 reveals that the empiricaldistributions of n A and n P are remarkably well fitted by a Beta distribution. This is exactly what one obtains inKirman and Föllmer’s ant recruitment model [
8, 9 ] , in which the Beta distribution emerges as the stationarydistribution describing a colony of ants preying on two distinct food sources. Such a distribution also emergesas the stationary distribution describing genetic populations between two competing alleles [
10, 11 ] . The keyingredient in these models is the competition between two different entities, be they food sources or geneticalleles. In Kirman and Föllmer’s ant model however, the two food sources are strictly equivalent and theresulting Beta distribution describing the fraction of ants at each source is necessarily symmetric, at odds withthe results obtained in the present setting. This motivates the asymmetric zones model introduced below.Another significant difference with Kirman and Föllmer’s original model is that the "food sources" here are notinexhaustible, as fish cannot reproduce at an infinite rate.These empirical results and observations motivate us to introduce a model extending Kirman and Föllmer’soriginal ant recruitment model to our context. In essence, one can think of the two cities as two distinct antcolonies that can obtain their food from either of the two zones. For each colony, the further fishing area isnecessarily less attractive, leading to the asymmetric character of the distribution. In addition, not being ina setting with unlimited resources, our model should take into account the fact that over-fishing may depletethe sea. III. A SIMPLE MODEL
Kirman and Föllmer’s original ant-recruitment model [ ] was successful at explaining a rather puzzling factwell known to entomologists [
12, 13 ] . Ants, faced with two identical and inexhaustible food sources tendto concentrate on one of them and occasionally switch to the other. Such intermittent herding behaviour isobserved in a variety of settings including choosing between equivalent restaurants [ ] , or financial mar-kets [ ] consistent with large endogenous fluctuations. In the model, at each time step a given ant mayeither (i) encounter another ant from the other inexhaustible food source and decide to switch to her peer’ssource (be recruited), or (ii) spontaneously decide to switch food sources without interacting. The drivingmechanism of the dynamics results from the trade-off between the intensity of the noise-term (cid:34) (spontaneousswitching), and that of the interaction term µ , see also [ ] .Here we present an extension of Kirman and Föllmer’s original model to account for exhaustible and asym-metric sources, notably aimed at accounting for some of the stylized facts presented in the previous section forfishing areas. Seeking to model fishermen exploiting a set of fishing areas, we imagine that boats follow thesame basic dynamics as the ants: if they initially fish within a certain zone, they may decide to move elsewhereeither because they see their peers fishing there, deciding to imitate them because they assume that their yieldis good, or spontaneously decide to move elsewhere randomly for the sake of exploration.As discussed above, our model has two major differences that depart from the original ant-recruitmentmodel. First, we consider that a fishing area has finite resources: fish reproduce until reaching a certain finitecapacity but they are also depleted by fishermen in the area (as in e.g. MacArthur’s models [
18, 19 ] ). As aconsequence, we assume that the random switching rate at which fishermen decide to depart from a givenarea depends on the fish population of that area. Note that this is very close in spirit to the modelling done inRef. [ ] , albeit that our model takes into account imitative behaviour in fishermen. The second difference withthe ant model is that we imagine two “colonies” instead of just one, corresponding to vessels based at the twodifferent fishing ports of Ancona and Pescara. Guided by the idea that fishermen prefer to go to areas close totheir own home port, notably to save fuel, we introduce an asymmetry between the fishing areas for each foodsource.The two ports, labelled A and P, have two distinct populations of fishermen, which may decide to exploit twofishing areas, S and S , with the fishermen from A preferring to fish at S and vice versa. One may of coursereasonably argue that this view is far too coarse-grained, and that there may be, for example, many differentfishing areas that are available close to each port. It is however possible to show under mild hypotheses thatthe two zones S and S in the model can be seen as the aggregation of a large number of smaller areas, withthe same dynamics, see Appendix B for details. For clarity, we shall define the model in discrete time, beforemoving into continuous time for analytical convenience.Without loss of generality, we focus only on the dynamics of fishing vessels at one of the two ports, sayAncona, as we assume that fishermen only interact with boats coming from the same city. We define now N A and N P as the number of boats based at Ancona and Pescara respectively, and let each of them decide to go toany of the two areas S and S . We denote m i ( t ) , with i =
1, 2, their respective fish populations at time t , andfurther assume that: • Boats only fish in one area each day (consistent with discussions with port authorities) and come backto that area if they don’t decide to switch to another one for the next day. • A vessel’s daily catch c i ( t ) is proportional to the amount of fish available in the area: c i ( t ) = β N A m i ( t ) with β/ N A ∈ [
0; 1 ] . • Fish reproduce at a multiplicative rate ν i , which we take to be equal to ν for both areas. • As a first approximation, fish do not travel from one area to the other. • The fish population within any area cannot exceed a carrying capacity K i , which is the maximal popula-tion that can be present within an area in the absence of fishing. This carrying capacity is the same forall areas, as we have taken all of them to be equivalent. Without loss of generality, we take K = K = Anecdotal evidence suggests indeed that the main interaction between people working in different boats happens at port in the fishingmarket or during informal conversation. Without changing our main conclusions, one could also allow for noise by drawing c i ( t ) from a given distribution centred about β m i ( t ) / N . This would allow us to introduce randomness into the fishing efficiency of each trawler, an interesting extension thatwe leave for further work. This constraint can be easily relaxed by e.g. adding a migration term where fish from 2 move to 1 at a certain rate and vice-versa. Inpractice, this would only tend to prevent the difference between the two fish populations from fluctuating too wildly.
Note that these definitions, which also amount to thinking of the fish population as consisting of the samespecies in both areas, are partially justified by our considering only trawlers, that therefore fish only veryspecific, shallow water dwelling species.We further define N A, i ( t ) the number of vessels from port A fishing at zone i at time t (and N P, i ( t ) respec-tively). The number of fishing vessels in each port is fixed, implying for all t : N A,1 ( t ) + N A,2 ( t ) = N A . Ourassumptions translate into following evolution for the fish population: m i ( t + ) − m i ( t ) = m i ( t ) (cid:2) ν g ( m i ( t )) − β (cid:0) N A, i ( t ) + N P, i ( t ) (cid:1)(cid:3) , (4)where the function g must satisfy g ( ) = g ( ) = ν m i g ( m i ( t )) models the amount of fish that are born between t and t + ν when the fish aren’t too plentiful, but goes to 0 as thezone gets saturated and cannot sustain any more fish. The simplest assumption one can make is that of logisticgrowth, leading to g ( m i ( t )) = − m i ( t ) . It follows that m i ( t ) ∈ [
0; 1 ] , ∀ t , where m = m = [ ] . The second term on the right handside corresponds simply to the decrease in fish population because of fishing activity.Furthermore, we assume that a fishing vessel based at A fishing at i can randomly decide to go elsewherewith probability (cid:34) A, i f ( m i ( t )) , where the function f satisfies f ( ) = f ( ) = + κ . Here, (cid:34) A, i controls thebase intensity of the noise, that can take a maximal value (cid:34) A, i ( + κ ) when the zone is depleted. Fishermenhave then a higher incentive to go elsewhere as their fishing yield decreases, and we highlight the preferenceof fishermen from A for zone 1 by setting (cid:34) : = (cid:34) A,1 = (cid:34) A,2 / C (cid:34) with C (cid:34) > S and S for a fisher from A. This allows us to have a higher spontaneous switchingrate S → S for fishermen from A.Besides this random switching rate, we add in the crucial element in our model, which is that agents imitateeach other. Each day, a fisherman randomly picks one of his peers and decides to imitate him with probability µ/ N , so that µ is the average fraction of boats deciding to take an imitation strategy at each step. In this case,the probability that a boat from A initially at zone S i decides to move to zone S j is given by: P A ( S i → S j ) = (cid:34) A, i f ( m i ( t )) + µ N A N A, j ( t ) N A − n A, i = N A, i / N A , we take the limit N A , N P → ∞ with N P / N A = C N fixed. Within this limit, wedenote n A = ( n A,1 , n A,2 ) and m = ( m , m ) and study the probability density ρ ( n A , n P , m ) of all the variables ofthe model. The equation one obtains is called the Fokker-Planck equation [ ] , also known as the Kolmogorovforward equation in applied mathematics and closely related to the Hamilton-Jacobi-Bellman equations de-scribing the optimal choice of a rational agent; it can be interpreted as the continuous limit of a Markoviantransition matrix. In the case of our model, the Fokker-Planck equation reads: ∂ t ρ = − (cid:34)∂ n A,1 (cid:2) C (cid:34) f ( m ) − n A,1 [ C (cid:34) f ( m ) + f ( m )] (cid:3) ρ + µ∂ n A,1 , n A,1 (cid:2) n A,1 ( − n A,1 ) (cid:3) ρ + (cid:2) ( n A,1 , m ) ↔ ( n P,2 , m ) (cid:3) − ∂ m (cid:2) ν ( − m ) − β (cid:0) n A,1 + C N ( − n P,2 ) (cid:1)(cid:3) m ρ + (cid:2) ( m , n A,1 , n P,2 ) ↔ ( m , 1 − n A,1 , 1 − n P,2 ) (cid:3) , (6)where the bracket [ x ↔ y ] is shorthand for the same expression where one replaces x by y .These equations fully close the model, which in our view represent the simplest setting for a system withlimited resources exploited by entities with a myopic exploration / imitation strategy. As they stand, however,they cannot be solved analytically. We shall now resort to a mean-field approximation to find a solution. As an interesting anecdote, we learned in [ ] that “Vito Volterra was born in the Jewish ghetto of Ancona in 1860, shortly before theunification of Italy, when the city still belonged to the Papal States”, and that “in 1925, at age 65, Volterra became interested in a studyby the zoologist Umberto D’Ancona, who would later become his son-in-law, on the proportion of cartilaginous fish (such as sharks andrays) landed in the fishery during the years 1905–1923 in three harbours of the Adriatic Sea: Trieste, Fiume and Venice. D’Ancona hadnoticed that the proportion of these fish had increased during the First World War, when the fishing effort had been reduced”. Thisled him to take interest in models that Alfred Lotka had first used to model very general population dynamics, and that we now apply,without knowing any of this at first, to the fish population dynamics at the ports of Ancona and Pescara. IV. MEAN-FIELD APPROXIMATION
Solving the Fokker-Planck equation (6) is a very difficult task, as the different terms that intervene take intoaccount the interactions between the proportion of fishermen in a given zone and the fish population in it. Ingeneral, these two quantities fluctuate in time, and the main difficulty lies in unravelling how these fluctuationsinteract. However, if one is interested in a very aggregated picture, one can simplify the problem significantlyby calculating the behaviour of the fish as if they were only subject to the averaged, fluctuation-free, actionof the fishing boats and vice-versa. This is what is known as the mean-field approximation , which allows toreplace, as a first step, the behaviour of the fish populations m i with their long-time averages. Taking Eq. (4)in the continuous time limit, the evolution of the fish population of, e.g., zone 1 follows:d m d t = m ( t ) (cid:0) ν ( − m ( t )) − β (cid:0) n A,1 + C N (cid:0) − n P,2 (cid:1)(cid:1)(cid:1) . (7)If we take the average of this equation, we expect the left hand side to be d (cid:69) [ m ] d t = This yields the following expression for the right hand side: (cid:69) [ m ] = (cid:149) − βν (cid:0) (cid:69) (cid:2) n A,1 (cid:3) + C N (cid:0) − (cid:69) (cid:2) n P,2 (cid:3)(cid:1)(cid:1)(cid:152) + , (8)where [ x ] + = x x > denotes the positive part of x . In particular, one can see that there exists a trivial extinctionline for the fish population for: ν = β (cid:2) (cid:69) (cid:2) n A,1 (cid:3) + C N (cid:0) − (cid:69) (cid:2) n P,2 (cid:3)(cid:1)(cid:3) , (9)which corresponds to the case where the reproductive rate of fish corresponds exactly to the rate at which theyare fished.We now do the same mean-field approximation the other way around to simplify the evolution of the fish-ermen. We insert Eq. (8) into the vessels’ dynamics by replacing the argument of f ( m i ) by the average, as f ( (cid:69) [ m i ]) : = f i , which amounts to saying that the boats only interact with the average behaviour of the fish.Choosing, to be precise, a linear function for f , i.e. f ( x ) = + κ ( − x ) , the average (cid:69) (cid:2) n A,1 (cid:3) can now be easilycomputed from Eq. (6) by setting the drift term to 0, as: (cid:69) (cid:2) n A,1 (cid:3) = C (cid:34) f C (cid:34) f + f = C (cid:34) ( + κ ( − (cid:69) [ m ])) + κ [ − (cid:69) [ m ] + C (cid:34) ( − (cid:69) [ m ])] . (10)Therefore, the mean-field approximation applied to the quantities m and n A have allowed to derive the twoEqs. (8) and (10), and therefore to obtain self-consistently the averages of these quantities. The next step isto obtain a fuller picture of the aggregate behaviour of the vessels within the mean-field approximation, andto obtain for example the probability density associated with it. A. Stationary solutions
Consistent with our mean-field approximation, we set m = (cid:69) [ m ] (resp. m = (cid:69) [ m ] ) in Eq. (6) to obtaina Fokker-Planck equation that describes the vessels. In the previous equation, vessels from the two portsinteracted indirectly through fishing in the same zone and depleting the fish in it, motivating all the boats in This average is taken with respect to infinitely many possible realizations of the stochastic process describing the evolution of our model.Nonetheless, because the model is ergodic we expect that this averaging is also true when one considers time-averages, meaning thatwe take a single realization of the process and look at the average of m in time. that zone to leave. Replacing the behaviour of the fish by its average amounts to neglecting this effect, and toan effective decoupling of the two variables n A and n P , as: ∂ t ρ = ∂ n A,1 J + ∂ n P,2 J , (11)where: J = − (cid:34) (cid:2) C (cid:34) f − n A,1 [ C (cid:34) f + f ] (cid:3) ρ + µ∂ n A,1 (cid:2) n A,1 ( − n A,1 ) (cid:3) ρ , (12)and where the transposition to find the definition of J is transparent. The two quantities J and J areprobability fluxes, and for example J (cid:0) n A,1 , t (cid:1) can be interpreted as the probability mass going from n A,1 + ∆ n to n A,1 for an infinitesimal ∆ n during an infinitesimal amount of time.The stationary state is found by setting J = J =
0, meaning that there is no probability flux in themodel, and solving the obtained equations for ρ . The decoupling of the two variables n A,1 and n P,2 allows oneto write the density as the product of two independent densities, as ρ ( n A,1 , n P,2 ) = ρ (cid:0) n A,1 (cid:1) ρ (cid:0) n P,2 (cid:1) , (13)with: ρ ( n A,1 ) = C n γ A,0 − (cid:0) − n A,1 (cid:1) γ A,1 − , ρ (cid:0) n P,2 (cid:1) = C n γ P,0 − (cid:0) − n P,2 (cid:1) γ P,1 − , (14)where C and C are normalisation constants and the γ parameters for the ρ distribution (the parameters for ρ can be easily deduced) read: γ A,0 = (cid:34)µ C (cid:34) f , γ A,1 = (cid:34)µ f . (15)Note also that full dynamical solutions ρ ( n A,1 , t ) , ρ ( n P,2 , t ) can be obtained in terms of hypergeometricfunctions, in the same spirit as [ ] , see Appendix A.For the reader unfamiliar with Fokker-Planck equations and stochastic processes, the following thoughtexperiment may help in understanding what we’ve stated mathematically above. The model in Sec. III can berun as a computer simulation where n A,1 and n P,2 correspond to numbers between 0 and 1 and can be knownat all times. It is then possible to run a very large number of simulations with the same initial conditions forthese two variables and to compute cross-sectional histograms for these variables at a given time-point, thatis the histograms of the same variable at the same time but through different simulations. The Fokker-Planckequation describes how these histograms will change in time, before eventually settling onto distributionsdescribed by Eq. (14). The model is however ergodic, and if we also take one single simulation and run it for avery long time, we can plot the histogram of, say, n A,1 ( { t i } ) at randomly sampled times t i and we will observea histogram described by the density ρ ( n A,1 ) .Our model thus replicates successfully the observed distributions shown in Fig. 2, and captures the qualitativebehaviour from Fig. 1. An example of numerical simulation of the model is provided in Fig. 3. For this Figure,we have set µ = κ = f , f ≈
1. We have then picked (cid:34) and C (cid:34) such as to match the values for γ and γ from Fig. 2. B. Dynamics and correlation functions
Within the mean-field model above, it is straightforward to show using tools from stochastic calculus (seethe appendices in Ref [ ] ) that the variable n A,1 follows a stochastic differential equation:d n A,1 d t = µ (cid:0) γ A,0 − (cid:0) γ A,0 + γ A,1 (cid:1) n A,1 (cid:1) + (cid:199) µ n A,1 (cid:0) − n A,1 (cid:1) η ( t ) , (16)0 n A ,1 − n P ,2 t m m n P > ( n ) Beta fit n A ,1 τ × − × − × − × − C − n A , ( τ ) Exponential fitSimulation
FIG. 3. Numerical simulation of our model. We have chosen the different parameters to obtain the same stationary Betadistribution as that observed in Fig. 2. Note also the similarity of the Figure on the left with plot (a) in Figure 1. We haveused the parameters (cid:34) = ν = β = C N = C (cid:34) = κ = n A,1 ( t ) and n P,2 ( t ) . The middle panel shows the fish populations m ( t ) and m ( t ) . The bottomleft panel shows the cumulative density function for n A,1 along with a Beta distribution fit. The bottom right panel shows theempirical correlation function as defined by Eq. (17) along with an exponential fit. Note that the fish populations oscillatearound the theoretical mean-field value m = with η a gaussian white noise of unit variance. This equation may be solved formally by integration, just as onewould for a standard ordinary differential equation. With this formal solution, it is then possible to computethe auto-correlation of 1 − n A,1 , defined as in Eq. (3), to find C − n A,1 ( τ ) = exp (cid:0) − µ (cid:0) γ A,0 + γ A,1 (cid:1) τ (cid:1) , (17)which is exactly what one sees from the data in Fig. 1 (d) and (e), provided one interprets the delta-peak at τ = n ( t ) = ( − σ ) n ( t ) + σξ ( t ) , (18)1 τ (days) − − − − C σ ( τ ) Ancona τ (days) − Pescara Direct exponential fitTheoretical decay
FIG. 4. Empirical correlation function C σ ( τ ) as defined in Eqs. (20) and (21). The solid black line is the theoreticalprediction given the estimations of γ and γ from the empirical probability distribution in Fig. 2 and from the subsequentestimation of µ using the exponential decay factor from Fig. 1. The reliable computation of σ depends of course on theproper estimation of these parameters, and we expect them to be noisy. Nonetheless, the agreement with theory, especiallyin the case of Pescara, is rather good. where n ( t ) is the “true” process and ξ ( t ) is a gaussian white noise of unit variance, then one can show directlythat the measured correlation function reads: C − ˜ n ( τ ) = δ ( τ ) + ( − σ ) σ C − n ( τ ) . (19)We have checked that the correlation function given in Eq. (17) agrees with our numerical simulations aswell; the results are shown in the bottom right panel in Fig. 3.The agreement is excellent both between the mean-field theory and the data, indicating that our modelcan correctly replicate the main dynamical features of real data from fishing dynamics. Furthermore, one candeduce the value of µ from the values of γ , γ and the decay factor in the exponential, that should match µ ( γ + γ ) . Using this formula, we find µ = · − for Ancona and µ = · − for Pescara.Note that the middle panel in Fig. 3 shows very interesting dynamics, with the fish population oscillatingabout its average value (cid:69) [ m ] = (cid:69) [ m ] = t ≈ m of fish close toAncona, and the population in that zone did not necessarily have the time to recover from that excess in thetime allowed by the simulation.More complicated correlation functions can also be computed, along the lines of [ ] , although they are moreprone to statistical noise. For example, using using the techniques described in detail in Ref. [ ] , one can showthat the polynomial defined by: σ A ( n A,1 ) = n − ( γ A,0 + ) γ A,0 + γ A,1 + n A,1 , (20)has an autocorrelation function that is exponential, meaning that C σ A ( τ ) = Cor (cid:0) σ A (cid:0) n A,1 ( t + τ ) (cid:1) , σ A (cid:0) n A,1 ( t ) (cid:1)(cid:1) verifies: C σ A ( τ ) = exp (cid:0) − µ ( + γ A,0 + γ A,0 ) τ (cid:1) , (21)and the same can of course be transposed to the variables indexed by P.2 β/ν E [ m ] Theory (mean-field)Simulation 0 1 β/ν "/µ γ Slope "κ Extinct phaseHigh noise phaseLow noise phase
FIG. 5. Left: Numercial results (same parameters as in Figs 3 and 4). The simulation was run for T = steps and with ν =
10. Note that the convergence of the simulation to the mean-field results from Eq. (23) improves as T or ν growlarger. Right: Phase diagram of the model. Note that the high / low noise frontier line could be extended in the extinctphase; indeed, without further ingredients, in the extinct phase boats move between the two areas as in the non-extinctphase even though they are not able to fish anymore. We have tested this prediction in Fig. 4. This correlator is necessarily more affected by measurement noise,because it is of order two in the n variables and because it depends on a reliable estimation of the γ and µ variables. Considering these limitations, the theoretical prediction is rather satisfactory when compared withthe data, especially in the case of Pescara. V. THE SYMMETRIC LIMIT
In general, the fixed point equations defined at the beginning of Section IV linking the averages (cid:69) [ m i ] withthe averages (cid:69) [ n i ] cannot be solved directly. Nonetheless, if one takes C N = f = f , withEq. (10) becoming: (cid:69) (cid:2) n A,1 (cid:3) = (cid:69) (cid:2) n P,2 (cid:3) = C (cid:34) C (cid:34) + (cid:69) [ m ] = (cid:69) [ m ] = (cid:26) − βν if β < ν β ≥ ν , (23)which has the intuitive interpretation that the population within a given area goes extinct if the fishing rate islarger than the reproduction rate of fish in the area. Figure 5 shows that the agreement of numerical simulationswith our mean-field analysis is excellent. One should note however that convergence may be slow when ν → f = f = + κβν . The parameters in Eq. (15) simplifyto yield: γ = ˜ (cid:34)µ C (cid:34) , γ = ˜ (cid:34)µ , (24)where we’ve dropped the A index as the parameters for both areas A and P are identical, and where we’ve set˜ (cid:34) : = (cid:34) (cid:128) + κβν (cid:138) .3 t n A ,1 − n P ,2 n P ( n ) Beta fit n A ,1 FIG. 6. Numerical simulations in the case γ >
1. The parameters are the same as that of Fig. 3, but with κ =
10 instead.The Beta fit is compatible with the predicted theoretical values from Eq. (24). Note that, in contrast with Fig. 2, we showthe probability density instead of the cumulative density function.
In this limit it is very clear that our mean-field model amounts to a modification of the original ant model [ ] , where the noise (cid:34) is augmented because of the sensitivity of the fishermen to the local fish populationby the factor given above, and where we have introduced an asymmetry between the two areas / food-sourcesthrough the parameter C (cid:34) .One would then typically expect to have always have γ > (cid:34) or κ are strong enough, one can have a crossover at γ = γ = n A,1 = n A,1 <
1, corresponding to γ > γ <
1, the empirical data shown in Figs 1 and 2 has γ < γ (cid:166) γ > ≈ − (cid:69) (cid:2) n A,0 (cid:3) of fishermen from Ancona fishing near Pescara. For the sake of completenessFig. 6 displays a simulation of this case, with γ well above 1. This last regime is qualitatively very differentto that with γ < VI. CONCLUSION
In this paper, we empirically analysed the distribution of the locations of fishing vessels in the two areas nearto Ancona and Pescara. By detecting to which port a vessel belongs, we computed the fraction of fishermenfishing in their home zone and looked at their statistical properties. We found that the empirical distributionfunctions are well approximated by asymmetric Beta distributions, and their auto-correlations by exponentials.Inspired by such evidence, we extended Kirman and Föllmer’s ants recruitment model to finite and asymmetricresources. We performed a numerical and theoretical analysis in the mean field approximation, and showedthat the auto-correlations and the stationary distribution of the fraction of fishermen appear to be respectivelyexponential and Beta distributed. Finally, we provided the phase diagram that separates a low and high herdingphase, as well as a fish extinction phase.We have further tested our dynamics by looking at higher order correlations that can be empirically com-puted. This signal appears to be very noisy and of low intensity but consistent with an exponential decay, with a4timescale compatible with that predicted by our model. Given the results that we have described, we are quiteconfident in our minimal model since it is able to reproduce surprisingly well the generic stylized facts withina limited though behaviourally sound framework. In particular, we have shown in Appendix B that while amulti-zone model (with more than two zones) would possibility be more realistic, the results for our two-zonemodel can be seen as the result of the aggregation of several zones, providing solid micro-foundations to ourapproach and justifying our looking at two aggregated zones for empirical analysis.
ACKNOWLEDGMENTS
We thank Jean-Philippe Bouchaud, Alexandre Darmon, Mauro Gallegati and Gianfranco Giulioni for fruitfuldiscussions and help in interpreting the data. This research was conducted within the
Econophysics & ComplexSystems Research Chair , under the aegis of the Fondation du Risque, the Fondation de l’Ecole polytechnique,the Ecole polytechnique and Capital Fund Management, and within the
New Approaches to Economic Challenges (NAEC) OECD programme on complex systems.5 [ ] Garrett Hardin. The tragedy of the commons.
Science , 162(3859):1243–1248, 1968. [ ] Stratis Gavaris. Use of a multiplicative model to estimate catch rate and effort from commercial data.
CanadianJournal of Fisheries and Aquatic Sciences , 37(12):2272–2275, 1980. [ ] Marianne Vignaux. Analysis of vessel movements and strategies using commercial catch and effort data from the newzealand hoki fishery.
Canadian Journal of Fisheries and Aquatic Sciences , 53(9):2126–2136, 1996. [ ] Peter M Allen and Jacqueline M McGlade. Dynamics of discovery and exploitation: the case of the scotian shelfgroundfish fisheries.
Canadian Journal of Fisheries and Aquatic Sciences , 43(6):1187–1200, 1986. [ ] Sorin Dascalu, Tudor Scurtu, Andreea Urzica, Mihai Trascau, and Adina Magda Florea. Using norm emergence inaddressing the tragedy of the commons. In
International conference on computational collective intelligence , pages165–174. Springer, 2013. [ ] Global Fishing Watch, Fishing Vessels Dataset, https://globalfishingwatch.org/datasets-and-code/vessel-identity , 2020. [ ] Mauro Gallegati, Gianfranco Giulioni, Alan Kirman, and Antonio Palestrini. What’s that got to do with the price offish? buyers behavior on the ancona fish market.
Journal of Economic Behavior & Organization , 80(1):20–33, 2011. [ ] Alan Kirman. Ants, rationality, and recruitment.
The Quarterly Journal of Economics , 108(1):137–156, 1993. [ ] José Moran, Antoine Fosset, Michael Benzaquen, and Jean-Philippe Bouchaud. Schrödinger’s ants: a continuousdescription of Kirman’s recruitment model.
Journal of Physics: Complexity , 1(3):035002, aug 2020. [ ] P. A. P. Moran. Random processes in genetics.
Mathematical Proceedings of the Cambridge Philosophical Society ,54(1):60–71, 1958. [ ] Sewall Wright. Statistical genetics and evolution.
Bull. Amer. Math. Soc. , 48(4):223–246, 04 1942. [ ] J. L. Deneubourg, S. Aron, S. Goss, and J. M. Pasteels. The self-organizing exploratory pattern of the argentine ant.
Journal of Insect Behavior , 3(2):159–168, March 1990. [ ] R. Beckers, J. L. Deneubourg, S. Goss, and J. M. Pasteels. Collective decision making through food recruitment.
Insectes Sociaux , 37(3):258–267, September 1990. [ ] Gary S. Becker. A note on restaurant pricing and other examples of social influences on price.
Journal of PoliticalEconomy , 99(5):1109–1116, 1991. [ ] David S. Scharfstein and Jeremy C. Stein. Herd behavior and investment.
The American Economic Review , 80(3):465–479, 1990. [ ] Robert J. Shiller and John Pound. Survey evidence on diffusion of interest and information among investors.
Journalof Economic Behavior & Organization , 12(1):47–66, August 1989. [ ] Thomas Lux. Herd behaviour, bubbles and crashes.
The Economic Journal , 105(431):881, July 1995. [ ] Robert MacArthur. Species packing and competitive equilibrium for many species.
Theoretical population biology ,1(1):1–11, 1970. [ ] Robert Mac Arthur. Species packing, and what competition minimizes.
Proceedings of the National Academy of Sciences ,64(4):1369–1371, 1969. [ ] Nicolas Bacaër.
Lotka, Volterra and the predator–prey system (1920–1926) , pages 71–76. Springer, 2011. [ ] Hannes Risken. Fokker-planck equation. In
The Fokker-Planck Equation , pages 63–95. Springer Berlin Heidelberg,1996. [ ] Milton Abramowitz and Irene A Stegun.
Handbook of mathematical functions with formulas, graphs, and mathematicaltables , volume 55. US Government printing office, 1948. Appendix A: Full dynamical solution
The goal of this section is to sketch a full dynamical solution for the dynamics of Eq. (16). We drop theindices A or P for clarity, obtaining:d n d t = µ ( γ − ( γ + γ ) n ) + (cid:198) µ n ( − n ) η ( t ) , (A1)a stochastic differential equation that corresponds to the following Fokker-Planck equation [ ] : ∂ t ρ = µ∂ nn ( n ( − n ) ρ ) − µ∂ n (( γ − ( γ + γ ) n ) ρ ) , (A2)with reflecting boundary conditions in n = n = [ ] , one can “diagonalize” this equation, writing it as: ∂ t ρ = (cid:65) ρ , (A3)with (cid:65) a Fokker-Planck operator that gives the right-hand side of Eq. (A2) when applied to ρ . It is in principlepossible to apply the same techniques as in Ref. [ ] to obtain a Schrödinger’s equation for an alternativefunction Ψ , that one could then use to compute ρ explicitly.However, one can also solve the eigenvalue problem (cid:65) ρ (cid:69) = (cid:69) ρ (cid:69) , so that the general solution for ρ reads: ρ ( n , t ) = (cid:88) (cid:69) λ (cid:69) ρ (cid:69) ( n ) e −(cid:69) t . (A4)In this setting, (cid:69) and ρ (cid:69) are respectively the eigenvalues and eigenvectors of the operator (cid:65) . These eigenvec-tors should also be normalized so that the integral of ρ is equal to 1.The problem therefore translates into finding functions ρ (cid:69) and numbers (or “energies”) (cid:69) that satisfy: µ (cid:69) ρ (cid:69) = µ∂ nn ( n ( − n ) ρ (cid:69) ) − µ∂ n (( γ − ( γ + γ ) n ) ρ (cid:69) ) (A5a) J (cid:69) ( ) = J (cid:69) ( ) = (cid:82) d n ρ (cid:69) ( n ) < ∞ , (A5c)with J (cid:69) ( n ) = µ∂ n ( n ( − n ) ρ (cid:69) ) − µ ( γ − ( γ + γ ) n ) ρ (cid:69) .In order to solve Eq. (A5a), we rewrite it as: n ( − n ) ρ (cid:48)(cid:48)(cid:69) + ( − γ − ( − γ − γ ) n ) ρ (cid:48)(cid:69) − ( + (cid:69) − γ − γ ) ρ (cid:69) = F ( a , b ; c ; n ) = (cid:88) k ( a ) k ( b ) k ( c ) k n k k ! , ( a ) k = k − (cid:89) i = ( a + i ) . (A7)Here the two linear independent solutions well defined around zero, see [ ] , are F ( a , b ; 2 − γ ; n ) and n γ − F ( a + γ − b + γ − γ ; n ) , where a , b are the solutions of: a + b = − γ − γ (A8a) ab = + (cid:69) − γ − γ . (A8b)Only the second solution cited above verifies the boundary condition at n =
0. Applying then an Eulertransformation on this solution leads to: ρ (cid:69) ( n ) = C (cid:69) n γ − ( − n ) γ − F ( − a , 1 − b ; γ ; n ) , (A9) The Euler transformation for the hypergeometric function states that F ( a , b ; c ; n ) = ( − n ) c − a − b F ( c − a , c − b ; c ; n ) . C (cid:69) a constant. This solution is well defined at n = (cid:90) d n ρ (cid:69) ( n ) = C (cid:69) (cid:88) k ( − a ) k ( − b ) k Γ ( γ + k ) Γ ( γ )( γ ) k Γ ( γ + γ + k ) k ! , (A10)with Γ the Gamma function. If 1 − a is a non-negative integer all terms in the series are non-zero. Usingthen ( x ) k ∝ Γ ( x + k ) together with the Stirling formula Γ ( x + ) ≈ (cid:112) π x x + / e − x for x (cid:29)
1, we find thatthe general term of the series converges to a constant when k → + ∞ and therefore that (cid:82) d n ρ (cid:69) ( n ) = + ∞ .Therefore, the condition that the functions ρ (cid:69) have a finite integral implies that there exists a positive integer k such that 1 − a = − k , and so also that b = − k − γ − γ and (cid:69) = − k ( γ + γ + k − ) . Since the numbers (cid:69) are discrete, and are indexed by k , we write ρ k : = ρ (cid:69) k .In conclusion, the eigenvectors ρ k and eigenvalues (cid:69) k are discrete and given by: (cid:69) k = − k ( γ + γ + k − ) (A11) ρ k ( n ) = C k n γ − ( − n ) γ − F ( − k , γ + γ + k − γ ; n ) , (A12)which allows then for a solution of the form given in Eq. (A4).There only remains to find the coefficients λ (cid:69) that depend on the initial condition. This can be done bytransforming the Fokker-Planck equation into a Schrödinger’s equation as in Ref. [ ] , noticing that the solutionsto said Schrödinger equation can be found in terms of the eigenvalues and eigenvectors ρ k , and one cantherefore find the coefficients λ (cid:69) by projecting the initial condition onto the orthogonal set of eigenvectors ofthe Schrödinger operators, see the Appendices in Ref. [ ] for a detailed technical explanation. Appendix B: A symmetric multizones extension
Here we present a very natural extension of our model to the general case of M symmetric zones with finiteresources. Without loss of generality we set C (cid:34) = n ( t ) = ( n ( t ) , . . . , n M ( t )) and m ( t ) = ( m ( t ) . . . , m M ( t )) ,where the index accounts for the zone, and call p j → i ( n ( t ) , m ( t )) the infinitesimal probability that an agentinitially present in zone j at time t moves to zone i at t + d t . It follows that the evolution of n and m is givenby: d m i ( t ) = m i ( t )( ν ( − m i ( t )) − β n i ( t )) d t (B1) p j → i ( n ( t ) , m ( t )) = n j ( t ) NM − (cid:2) (cid:34) f (cid:0) m j ( t ) (cid:1)(cid:3) + µ N n i ( t ) n j ( t ) . (B2)Introducing for simplicity h ( n , m ) = m ( ν ( − m ) − β n ) , the joint density ρ of the variables ( n ( t ) , m ( t )) evolvesaccording to the following Fokker-Planck equation: ∂ t ρ ( n , m ) = − (cid:88) i ∂ m i ( h ( n i , m i ) ρ ) + (cid:88) i (cid:54) = j (cid:128) ∂ n i − ∂ n j (cid:138) (cid:130) (cid:34) f (cid:0) m j (cid:1) M − n j ρ (cid:140) + µ (cid:88) i (cid:54) = j (cid:128) ∂ n i n i − ∂ n j n i (cid:138) (cid:0) n i n j ρ (cid:1) . (B3) As the hypergeometric function is symmetric with respect to its two first arguments, we restrict our analysis to the first one only. (cid:69) [ n i ] = M , (cid:69) [ m i ] = (cid:129) − β M ν (cid:139) + (B4)with ( x ) + the positive part of x . This again shows the existence of an extinction regime whenever β = M ν . Inwhat follows we assume that β/ν < / M , to study the behaviour of the system outside of extinction.When κ = n is a Dirichlet distribution with all parameters equal to ( (cid:34)/δ ) , namely: ρ n ( n , . . . , n M ) = (cid:130) M (cid:89) i = n (cid:34)/µ i (cid:140) { (cid:80) Mi = n i = } . (B5)As argued previously, whenever the noise-level is coupled to the fish population with κ >
0, we postulatethat the solution can be approximated by a Dirichlet distribution with all parameters set to ( ˜ (cid:34)/µ ) with ˜ (cid:34) = f ( − β M ν ) (cid:34) .The Dirichlet distribution has one key property: (cid:80) i ≥ k n i follows a Beta distribution with parameters ( k ˜ (cid:34)/µ , ( M − k ) ˜ (cid:34)/µ ) , corresponding to the stationary state of our two-zone model. We have also checkedthat the mean-field approximation of Eq. B1 follows the same type of property: the variable (cid:80) i ≥ k n i ( t ))