Author’s Accepted Manuscript Models for Allee effect based on physical principles Renato Vieira dos Santos, Fabiano L. Ribeiro, Alexandre Souto Martinez www.elsevier.com/locate/yjtbi
PII: DOI: Reference:
S0022-5193(15)00411-7 http://dx.doi.org/10.1016/j.jtbi.2015.08.018 YJTBI8325
To appear in: Journal of Theoretical Biology Received date: 13 March 2015 Revised date: 21 August 2015 Accepted date: 24 August 2015 Cite this article as: Renato Vieira dos Santos, Fabiano L. Ribeiro and Alexandre Souto Martinez, Models for Allee effect based on physical principles, Journal of Theoretical Biology, http://dx.doi.org/10.1016/j.jtbi.2015.08.018 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Models for Allee effect based on physical principles Renato Vieira dos Santos1 , Fabiano L. Ribeiro2 and Alexandre Souto Martinez3 1
UFLA - Universidade Federal de Lavras, DFI - Departamento de F´ısica, CEP: 37200-000, Lavras, Minas Gerais, Brazil.,
[email protected] 2 UFLA - Universidade Federal de Lavras DFI - Departamento de F´ısica, CEP: 37200-000, Lavras, Minas Gerais, Brazil.,
[email protected] 3 USP - Universidade de S˜ao Paulo, FFCLRP - Faculdade de Filosofia, Ciˆencias e Letras de Ribeir˜ao Preto, Departamento de F´ısica e Matem´atica, Av Bandeirantes 3900, CEP: 14040-901, Ribeir˜ao Preto, S˜ao Paulo, Brazil.,
[email protected] August 28, 2015 Abstract We propose some models of single species with Allee effect based on physical principles. A method is used to obtain the expression for the per capita growth rate (a macroscopic information) starting from the characteristics of interactions between the individuals (a microscopic information). We assume that the agents in a model of a single species interact according to the distance between them. Moreover these agents must (i) cooperate with their nearest neighbours, (ii) compete with neighbors at an intermediate distance, and (iii) being indifferent to those who are far away. Using these assumptions and based on fundamental physical principles, we find what appears to be a new way of establishing models of single species with Allee effect.
1
Introduction
species with very small populations did not compete with each other. He noted, for example, that aggregation has a positive effect on survival of animals of these species. When isolated from their community in extinction, these animals died more quickly than their fellow species.
“An instinct is a more or less complicated activity of an organism which is acting (1) as a whole rather than as a part; (2) as a representative of a species rather than as an individual; (3) without previous experience; and (4) with an end or purpose of which it has no knowledge.” [1]
Allee demonstrated this effect in his book with real aggregation experiments, contradicting the paradigm of the logistic growth model. In the logistic model, a population has a slow growth in the beginning, a phase of rapid acceleration and finally a phase of stationary growth. Allee observed that in a small population such growth may exist but there may be also the extinction of the species. Everything depends on a critical point that is mathematically determined as the point of system instability. Below this critical point - the so-called Allee
These were the words of the great biologist and ecologist Warder Clyde Allee in the first chapter of his book “Animal Aggregations”, when he began to talk about instinct [1]. Thanks to the studies of Allee, a pattern observed in some communities of organisms received its name, Allee Effect. Allee, being an observer of animal behavior, realized that individuals of some 1
threshold - the population vanishes and above it the population has a positive growth phase. In other words, a population with the Allee effect may have a threshold population size below which the population goes extinct deterministically. Basically, the “Allee effect” has two manifestations [2]:
All mentioned models are described as dN = N G(N ), dt
(3)
where G(N ), the per capita growth rate function, is λ λe−αN (t) −1 for the Ricker model and (1+αN −1 (t))β for the Hassel model. None of the mentioned mod1. There is a positive association between some els has Allee effect because G(N ) are decreasing component of fitness and population size (e.g., functions of N. When a population faces the Allee increasing part followed by a fecundity, viability of the existence of groups, effect, G(N ) has an 2 G(N ) etc). > 0) at low N (see figdecreasing part ( d dN 2 dG(N ) 2. There is a positive association between per ure (1)). If G(N = 0) < 0 and dN N =0 > 0, capita population growth and population size. we have strong Allee effect. If If G(N = 0) > 0 dG(N ) Allee observed this association in colonies of bac- and dN N =0 > 0, we have weak Allee effect [2]. teria, fungi and endangered species such as turtles Many expressions have been used in the literature or birds [2]. These species also realize that from to describe Allee effect, most of them proposed phethe same point on may be time to start compet- nomenologically, of which the following two are the ing to preserve the same small size. The growth most frequently used: of the species may not be interesting because there N a N may be the attraction of predators. Few in num− G(N ) = k 1 − , (4) K K K ber, these species will be less apparent and survive longer. Allee effects can also occur when populations perform poorly when density is low because and of inbreeding, low probability of finding mates, cola+τ N lapse of social systems, lowered capability to deter 1− . (5) G(N ) = k 1 − K N +τ predators, etc. Some recent experimental results of population which exhibit Allee effect are in [3, 4, 5]. In both models, k is the maximum per capita Several models based on different assumptions growth rate in absence of the Allee effect (the inhave been proposed to describe the population trinsic growth rate), K is the carrying capacity, a is dynamics. These models are known by the the Allee threshold (see figure (1)), and τ in Eq. (5) names of their authors: Ricker, Gompertz, Hassell, is a non-negative parameter affecting the overall Beverton-Holt, Maynard-Smith, Malthus, Verhulst shape of G(N ). and so on [6]. For instance, the continuous version of the Ricker model [7], dN (t) = N (t) λe−αN (t) − 1 , (1) dt where N (t) is the population size at time t, λ is the intrinsic growth rate and α is the per capita competitive effect, has been used to model successfully animal populations growth [8]. Another example, the Hassel model [9] λ dN (t) = N (t) − 1 , (2) Figure 1: Per capita growth function G(N ): withdt (1 + αN (t))β out Allee effect (black dotted line), with weak Allee provides the best fit to plant data. Frequently β effect (blue dashed line), and with strong Allee efis equal to one, so the Hassell model becomes the fect (red dotted line). Allee threshold is the black Beverton-Holt model in most plant species [10]. dot to the left. 2
In this paper, motivated by physical concepts, we propose a way to obtain G(N ) starting from assumptions regarding the interactions between the agents of the species which are spatially distributed. Specifically, we assume that the interaction between the individuals depends only on their distance. That is the case, for instance, of the yeasts described in [11]. The cooperation between these yeasts decays with the distance because the common good produced by a single cell diffuse away to be consumed by other cells in the population. We will consider that the strength of the interaction between two individuals is described by a function f (r), where r denotes the distance between them. Moreover, the signal of f (r) gives the kind of interaction between two individuals: if positive, they cooperate; if negative, they compete. The function f (r) has information about the microscopic interactions, given that it describes the way that the individuals interact. In this present work will be presented, as one of the results, the correspondence
Figure 2: Required general form of strength of interaction between two individuals f (r). Cooperation exists for r < rc and competition for r > rc . For r = rm competition is maximum and f (r) → 0 when r → ∞.
2
G(N ) from physical principles
iii) For r > rc , there is a rm > rc such that competition is maximum and f (r) → 0 when r → ∞.
Consider a population formed by N (t) individuals at time t. This population is embedded by a hypervolume VD in a D−dimensional euclidean space. Consider that the interaction strength between any two individuals, say i and j, depends only on the distance rij = |rij | between them, where rij ∈ RD . More specifically, we will consider that the interaction strength between a pair of individuals will be described by the microscopic model f (rij ). The γ , where γ is the decay particular case f (rij ) = 1/rij parameter, was studied in [12]. Other specific cases of this microscopic model was studied in [13]. Following the ideas introduced by Mombach et al [12] and reworked by Ribeiro [14, 15], the replication rate Ri of the i−th agent of the population is composed by two factors: i) the self stimulus to replicate, and ii) the stimulus - the field - from the other individuals of the population. In this way, the replication rate of the agent i can be written as f (rij ). (7) Ri = k +
Commonly used models does not exhibit the property iii, which is not physically plausible. We propose and study some models that exhibit all the properties above. The general form of f (r) must be as shown in Figure (2).
The parameter k is the self-stimulus of the individuals; it is identical for all agents and does not depend on the population size. This parameter represents all the replication aspects of the individual which do not depend on the other individuals of
G(N ) ⇔ f (r).
(6)
That is, starting from G(N ) (the macroscopic information) we obtain f (r) (the microscopic information) and vice versa. We will use this correspondence to obtain f (r) for various models (various G) established in the literature. To our surprise, to the best of our knowledge, none of the models commonly used to describe the Allee effect has certain properties considered desirable for the function f (r). Let us define rc as the critical distance. The desirable properties for a model that describes the Allee effect are: i) For r < rc , there must be cooperation. ii) For r > rc , there must be competition.
j=i
3
the same. That is, the population is homogeneously distributed when DF = D. As f (r) is a function that depends only on the distance r, it is interesting to use hyper-spherical coordinates to solve the integral in Eq. (8). Identifying dD ri = rD−1 dΩD dr, assuming periodic boundary conditions, and using (9), Eq. (8) becomes
Rmax f (r)rDF −1 dr. (10) G(N ) = k + ρ0 ΩD
the population. In fact, k is the growth rate of Malthus model, as we present in the following (see subsection (4.1)). Given ρ(r), the density of individuals at the position r, and ri , the position of the agent i, we can show (see appendix A) that the per capita growth rate G(N ) = (dN/dt)/N can be computed from the microscopic model by G(N ) = k +
N
1 dD rρ(r)f (|r − ri |), N i=1 VD
(8)
Here ΩD = dΩD represents the solid angle, which values are: Ω1 = 2, Ω2 = 2π and Ω3 = 4π. This radial symmetry assumption implies some limitations. In this case, only individuals inside the cluster of individuals are modeled. Any edge event is neglected as a consequence of the exclusive radial dependence of f. This radial dependence comes from the corpuscular character of the interacting agents. Since agents are treated as particles, there is not much sense in imagining that the interaction between them is anything but a radial interaction. The extremes of the integral, Rmin and Rmax , are the minimum and maximum euclidean distance between two cells. The distance Rmin can be, for instance, the diameter of the cell, and the distance Rmax is the necessary radius to embedded all the cells of the population. The dependence on N in Eq. (10) comes from Rmax because, as was showed 1 in [12, 15], Rmax = [DF N/(ΩD ρ0 )] DF . Extensive analysis of the population behavior given different fractal dimensions were studied in [15]. In conclusion, Eq. (10) represents a systematic way to obtain the per capita growth rate of the population with fractal distribution from the microscopic model - the interaction function - f (r). In this way, we obtain from Eq. (10) the population properties by a non-phenomenological approach, in opposition to what is usually done in the context of population modeling.
allowing us to obtain information about the macroscopic behavior (per capita growth rate G(N )) from the microscopic information f (r). In short, Eq. (8) shows that it is possible to get the macroscopic growth behavior of the population if we have: a) the spatial distribution of the individuals, characterized by ρ(r); and b) the way that the individuals interact, characterized by f (|rij |).
2.1
Population structure with fractal dimension
There are many examples of microorganisms which its populations exhibit fractal dimension, as is the case of tumors [16, 17, 18, 19] and colonies of bacteria [20, 21]. Populations which exhibit fractal-like structure, are self-similar. That is, its parts show the same statistical properties at many scales. In the case of tumors, the fractal structure of the population is formed due to the irregular reproduction patterns of the cells. In addition, the fractal dimension of the cancer is strongly related to its aggressiveness [22]. In the case of colonies of bacteria, the fractal-like structure is due, among others, to the step by step repeated formation of micro-bursts [20]. Let us assume that the population forms a structure with fractal dimension DF . Then the structure formed by this population must have selfsimilarities and the density of agents can be written in terms of a generic radius r. In this case, the density of individuals is (see [23]) ρ(r) = ρ(r) = ρ0 rDF −D ,
Rmin
3
The G(N ) ⇔ f (r) correspondence
(9) When dealing with populations growths, usually the only quantities that can be measured are the macroscopic properties of the population. For instance, it is relatively easy to obtain at the mi-
where ρ0 is a constant. The constant ρ0 represents the density of cells when the dimensions of the population structure and the hyper-volume are 4
4.1
crobiology laboratory some macroscopic quantities as: total number of individuals as a function of time; the density of individuals; and the per capita growth rate (which can be obtained fitting the population growth data curve). However, the microscopic properties of the system, for instance, the way that the cells interact, is very difficult to access experimentally. The method that we propose gives a route to get the microscopic properties (hard to get) from the macroscopic information (usually available experimentally). It is possible solving Eq. (10) in order to isolate f (r). As shown in details in appendix B, dG(N ) , (11) f (r) = dN D
Malthus model
The simplest population growth model is the exponential growth, described by the Malthus model. In this model, each individual reproduce at a fixed rate with G(N ) = constant = k and, from Eq. (11), f (r) = 0. Individuals of the population did not interact. From Eq. (3) one gets dN/dt = kN. In according to what was mentioned earlier, the self stimulus of the cell, k, is the Malthus constant growth rate.
4.2
Verhulst model and generalizations
The Verhulst model has a saturation regime in the population growth. The model can be represented by the per-capita growth rate G(N ) = k(1 − N/K), N =μr F where K is the carrying capacity of the population. ΩD with μ ≡ D ρ0 . The quantity G(N ) is easy to Using Eq. (11) one gets F access experimentally and so, from it, one can calk (12) f (r) = − = −constant. culate, by Eq. (11), the interaction distance depenK dent function between any two individuals in the Verhulst model can be interpreted as a kind of mean population. field model where the interaction between the indiImportantly, several functions G(N ) lead to the viduals happens regardless of the distance between same function f (r). It is sufficient to change k. This them. Moreover, given k > 0, f (r) is always negameans that for a significant change in the radial tive, which implies that the individuals are always function f (r), it is necessary that there are non in competition. linearities in dN/dt = N G(N ). This makes sense The generalized model is defined as G(N ) = from a physical point of view because the interac- k (1 − (N/K)q˜), which is also known by Richards q˜ tions (non linearities) are responsible for affecting model [24, 25]. Here, q˜ is the generalization paramthe function f, by definition. eter: for q˜ → 0 the model retrieves the Gompertz It is worth noting that the requirement of unique- model; for q˜ = 1, the model retrieves the Verhulst ness among G and f and finiteness of the popu- model. For more details see [26]. Using Eq. (11), lation for long times imply a minimum model for one gets an algebric power dependence the Allee effect. Under these conditions, we would 1 have G(N ) = kN − βN 2 , where the first term is (13) f (r) ∝ − D (1−˜q) . F r the element of lowest order in N to which there is This is exactly the form of the interaction uniqueness between G and f, and the second term strength chosen by Mombach et al which give rise ensures finite population in the long term. to the Richards model in [12]. Once more, f (r) is always negative, and individuals are always in 4 Applications of the method competition. We test a wide gamma of possibilities with a generalized Verhulst model proposed in [27]: In this section, we present the applications of β γ Eq. (11) to some well known models. We start N α−1 . 1− G(N ) = kN with Malthus and Verhulst models and then we K analyze their generalization, the Richards model. Finally we study continuous versions of the Ricker Here, α, β and γ are positive parameters. There is only competition again. and Hassel models. 5
4.3
Ricker model
intense competition neglecting distant individuals hinders the formation of groups. In the introduction we mentioned the Ricker model The Hassel model inspired by Eq. (54) of ap(1). To obtain this model from first principles, it is pendix C returns the function (17) with n renorsufficient to assume that individuals are distributed malized by n/s, as in the Ricker model and the in n sites of a lattice according to Poisson and that interpretation is analogous. in each site containing k = 1 individuals, b offspring arise (see appendix C). Parameterizing the model according to [8]: 5 Common models with Allee N dN = bN e− n − N. dt
effect
(14)
In this section we study some classical models that exhibit Allee effect for future comparisons with the models that will be proposed.
N
Placing G(N ) = be− n − 1 in Eq. (11), we have b μ DF f (r) = − e− n r n
(15)
5.1
Models with multiplicative Allee
where f (r) < 0 for all r > 0, resulting only in comeffect petition. Rising b and n, the strength of competition (given by |f (r)|) and the range of interactions We consider the following family of models given increase, respectively. These facts are consistent by differential equations of the general form with the interpretations given in the appendix C. N If we consider the continuous Ricker model inG(N ) = k 1 − A(N ). (18) K spired by Eq. (53) in appendix C, we get Eq. (15) with n renormalized to n/(1 − C) ≡ n/s, 0 ≤ s ≤ 1 The function A(N ) incorporates an Allee effect. denoting a microscopic parameter related to the inSince A(N ) multiplies the logistic term, we have a tensity of struggle between individuals. Decrease of multiplicative Allee effect [6, 28], in contrast to the the intensity of the dispute is equivalent to increasadditive case found in [29, 30, 28]. ing the amount of sites n.
4.4
5.1.1
Hassel model
The simplest model
A simple phenomenological model is the one with The Hassel model can be obtained considering that per capita growth rate given by individuals are initially distributed in accordance with the negative binomial distribution, as shown N N −1 , (19) G(N ) = k 1 − in the appendix C. From Eq. (50), in the continuous K K case, bλλ+1 with A(N ) = N/K − 1. The per capita growth G(N ) = −1 (16) λ+1 rate is negative when N < K , which implies that λ+ N n the population extinguish when that is the case. and, from Eq. (11), When N > K , the population presents a positive growth rate and G increases as the population sizes b(λ + 1)λλ+1 increases. In a specific size, G became to decrease f (r) = − . (17) λ+2 and then vanishes at N = K. That is the stable n λ + nμ rDF equilibrium point. This function supports only competition. Following Using Eq. (11), one gets interpretation of [8], λ is associated with the degree of clustering of the individuals. Large values of λ (20) f (r) = −arDF + b, imply less grouped individuals. Increase of λ in 2kμ K +K > 0. The Eq. (17) has two effects: increase the strength and where a ≡ KK > 0 and b ≡ k K K decrease the range of the interaction. We see that interaction function f (r) can be expressed by two 6
terms. The first one (on the right hand side of 5.2 Model with additive Allee effect the Eq. 20)) is negative and then is responsible for The following model with additive Allee effect will the competition. The second term is positive and be considered [29]: responsible for the cooperation. While the cooper ative term is constant, the competitive term grows N km G(N ) = k 1 − . (22) − DF as a power r . The result is a positive interaction K N +τ (cooperation) for all individuals which is stronger when individuals Above the critical dis- Eq. (22) can be rewritten as 1 are closer. K +K DF , we have competition, as tance rc ≡ kN N 2μ G(N ) = 1− (N − N− ) (23) N +τ N+ shown in Fig. (3). As limr→∞ f (r) = −∞, then Allee property iii mentioned in the introduction is not Allee satisfied. with N± = 12 K − τ ± k1 [k(K + τ )2 − 4mK] . This model have been used to describe multiple Allee effects [31, 29, 32]. In the case described by (23), a double Allee effect. From Eq. (11), km k k f (r) = 2 − K with limr→∞ f (r) = − K , (μrDF +τ ) violating property iii.
6
Proposed models
In this section we develop two simple models that satisfy the three properties formulated in the introduction. In addition to satisfying the physical requirements for f (r), they are motivated by first principles summarized in the appendix C. From these principles, by demanding that only in sites containing k = 2 individuals may be reproduction and using a) Poisson and b) negative binomial distributions, we obtain models similar to the models of Ricker and Hassel, respectively, contemplating Allee effect. Using Eq. (47) with φ(k) = bδk,2 and
Figure 3: The plot of f (r) for the model (19) which displays Allee effect. Note two regimes depending on the distance: Cooperation at distance r < rc and competition for r > rc .
pk = (Nt /n)k e(−Nt /n) /k!, 5.1.2
Another model
we obtain the following power -Ricker model:
Another widely used model [2] is described by Eq. (5), which is Eq. (18) with A(N ) = 1 − Na+τ +τ where a is the Allee threshold and τ affects the with overall shape of G(N ). From Eq. (11):
f (r) =
bN 2 − N dN = e n −N dt 2n
(25)
bN − N e n − 1. (26) 2n For our purposes, the constant term on the righthand side of Eq. (26) can be generalized to an arbitrary value. Thus, we consider the following subtle generalization:
k a(K + τ ) − μrDF μrDF + 2τ + Kτ K (μrDF + τ )
(24)
2
(21) k . Property iii is violated with limr→∞ f (r) = − K again.
G(N ) =
G(N ) = 7
bN − N e n − g. 2n
(27)
For Allee effect with an equilibrium solution other that N = 0, we must have b > 2ge. For g = 0 we have weak Allee effect. Corresponding to Eq. (27), we get b − μ r DF n − μrDF e n 2n2
(28) Interaction Intensity f (r)
f (r) =
From Eqs. (31) and (32) we see that the effect of decreasing λ is to increase the intensity of cooperation at the cost of reducing its range. This fact is consistent with the interpretation given in the appendix C1 and it is well represented in Fig. (5).
with b f (0) = 2n
and
D1 F n rc = . μ
(29)
Fig. (4) illustrates the dependence of f (r) with b. The figure is similar for power Hassel model defined below.
0.6 0.3 0 0
0.6 Interaction Intensity f (r)
0.9
1
2
3
4
5
6
Distance (r)
0.5 0.4
Figure 5: Interaction strength for power -Hassel model for λ = {1,2,10} with {black,blue,red} curves, respectively. Other parameters: μ = 1/DF , n = 10, b = 10, DF = 3.
0.3 0.2 0.1
Fig. (6) shows the dependence of f (r) with DF for the power -Ricker model. Similar behavior occurs for power -Hassel. 0 1 2 3 4 5 6 An efficient way to investigate the balance beDistance (r) tween cooperation and competition is by the calculation of the Social Cooperativity (SC), defined Figure 4: Interaction strength for power -Ricker by
∞ model for b = {1,5,12} with {black,blue,red} f (r)dr. (33) SC ≡ curves, respectively. Other parameters: μ = 1/DF , 0 n = 10, DF = 3. If SC > 0 (SC < 0) we have more cooperation (competition) than competition (cooperation). Similarly, for the negative binomial distribution, If SC = 0, the community cooperates as much we have the power -Hassel model: that competes. Satisfy property iii means higher chances of having a finite SC. For the power -Ricker b λλ+1 (λ + 1)N (pRi) and power -Hassel (pHa) models, SC will be, − g (30) G(N ) = 2n λ + N λ+2 respectively n b(DF − 1)Γ D1F with , (34) SCpRi = 1 2DF2 n nμ DF bλλ+1 (λ + 1)n λn − (λ + 1)μrDF (31) f (r) = λ 3 2 (μrDF + λn) nμ rDF + λ b(DF − 1)Γ D1F Γ λ − D1F + 2 , SCpHa = μ D1 and F 2DF2 λ2 nΓ(λ) λn D1 (35) F n λ b λ+1 f (0) = , rc = . (32) 1 Low λ is related to clustering. 2n λ μλ+1 0
8
where Γ(x) is the gamma function. It is clear that large values of b (number of offspring per site) favors and large values of n (number of sites) disfavors cooperation. The ratio of equations (35) and (34) is given by Γ λ − D1F + 2 SCpHa (36) = 1 −2 SCpRi 1 DF Γ(λ) λ
result: n b
Interaction Intensity f (r)
0.75 0.5 0.25 0 5
7.5
t-Stat 6.17 8.77
P-Value 6.25 × 10−6 4.20 × 10−8
0.781 −50.0 −46.87 0.802
In the tables above we see the estimated values of the parameters b and n and some goodness of fit measures such as R Squared, AIC (Akaike information criterion [34]) and BIC (Bayesian information criterion [34]). These two latter are characterized by penalizing the number of parameters of a model. More negative values indicate better models. In light of the biological interpretation, it is interesting to note how close to 1 is the parameter b, with little standard error. This is consistent with the fact that the musk oxen have one baby per birth. For a biological interpretation of n in the continuous time context, we make an expansion in Taylor a in the variable ≡ 1/n 1, giving, for series of g=0: b G(N ) ≈ N (1 − N ). 2
1
2.5
Std Error 32.71 0.113
AdjustedRSquared AIC BIC RSquared
and depends only on the parameters DF and λ. For DF > 1, SCpHa /SCpRi > 1 for all λ with limλ→∞ SCpHa /SCpRi = 1. Furthermore, the value of λ to which Eq. (36) is maximum is λ∗ = 1 and for this λ, Eq. (36) gives Γ(3 − 1/DF ). More crowding (λ → 1), more cooperation (maximum of (36)).
0
Estimate 201.872 0.98891
10
Distance (r)
Solving G(N ) = 0 results in N ∗ = 1/ = n, which Figure 6: Interaction strength for power -Hassel is the carrying capacity in this approximation. In model for DF = {1.5,2,3} with {black,blue,red} discrete time n is the number of sites and in concurves, respectively. Other parameters: λ = 1, tinuous time n is related to the carrying capacity. The data and the curve fit are shown in Fig. (7). n = 10, b = 10 and μ = 1/D . F
7
Application to data
In this section we apply the various models presented in table below to musk ox data directly obtained from [33]. Model G(N) = N bN a power -Ricker 2n exp − n − g b N b Ricker 2n exp − n − g c Basic Allee k (1 − QN ) (Q N − 1) λ+1 (λ+1) b λ d Hassel 2n (λ+ N )λ+2 − g n λ+1 (λ+1)N b λ e power -Hassel 2n (λ+ N )λ+2 − g n
Figure 7: Musk ox data fitted by power -Ricker model with g = 0.
a with g = 0 to the data Adjusting the function using Mathematica software, we have the following
For the power -Ricker model with g = 0 there are 9
no finite equilibrium solution and G(0) = 0 (No ure (8) shows f (r) for the two estimated models. strong Allee effect). Fitting the model with g = 0 we obtain: b n g
Estimate Std Error 151.085 51.98 0.826 0.181 −0.000124 0.00012
t-Stat P-Value 2.91 0.00941 4.55 0.000247 −1.04 0.312
AdjustedRSquared 0.777 AIC −48.783 BIC −44.605 RSquared 0.809 Observing the P-Value of g, there is no reason to consider it different from 0. This implies that there are no evidence for strong Allee effect for the data analyzed. Figure 8: Estimated f (r) for power-Ricker (pRi) e with λ = 1, 2, · · · , 6. The and power-Hassel (pHa) models. We fit the function best result was for λ = 1:
b n
Estimate 1.13 354.708
Std Error 0.126 81.713
t-Stat 8.990 4.341
P-Value 2.84 × 10−8 0.000352
AdjustedRSquared 0.79 AIC −50.83 BIC −47.69 RSquared 0.81 This was the best result found considering all the analyzed models. The fact that λ = 1 being the best value to fit the power -Hassel model is consistent with the fact that muskoxen live in herds [35]. As discussed in the appendix C, small values for λ involve the existence of agglomerates. e in Taylor series in the variable ≡ Expanding 1/n 1 as before, the carrying capacity for the power-Hassel model will be N ∗ = λn. Therefore, for musk ox data, N ∗ ≈ 350. It is known that a musk ox population fluctuating between 200 and 300 animals has occupied the northeast coast of Devon Island for several hundred years [36, 37]. If a condition to have good mathematical models is the explicit manifestation of the underlying reality being modeled, it seems that we are on track. With the estimated parameters and some more information we can estimate f (r) for the two models, power-Ricker (pRi) and power-Hassel (pHa). Considering Df = 2, n = 250 and that the area of the Devon island is approximately 50 km2 , we have ρ0 ≈ 250/50 individuals per km2 and μ ≈ 5π. Fig-
We see that the distance from which the competition start is approximately R ≈ 3.5 km for the two models. This distance corresponds to an area of A = πR2 ≈ 40 km2 . Considering the density ρ0 of 5 individuals per km2 , we have that approximately 200 muskox can live cooperatively. Typical results for the other models are well represented by the example below. c we have: For function k Q Q
Estimate −0.097 0.0013 −0.0055
Std Error 0.0312 0.000274 0.0046
AdjustedRSquared AIC BIC RSquared
t-Stat P-Value −3.116 0.006 4.788 0.000147 −1.20 0.244
0.78 −49.18 −45.00 0.81
AIC and BIC, by penalizing the number of parameters, indicate the inferiority of this curve fitting (three parameters) to the previous (two parameters). Far more significant, there is no reason to believe that Q is not 0 because your PValue(=0.244) is too large.
10
8
Conclusion
In the present work it was developed a framework which can be a applied to any kind of population
in which the intensity of the interaction (cooperative or competitive) between the individuals depends only on the euclidean distance between the individuals. Based on this approach we propose two models that incorporates Allee effect based on physical first principles. Considering that individuals interact with intensity f (r) depending on the distance r between them, we find a connection between this function (a microscopic information) and the per capita growth rate G(N ) (a macroscopic information). We tested the results in several established models for the description of the Allee effect and found that none of them satisfies certain properties considered desirable for f (r). It is interesting that a reasonable assumption on f (r) is not consistent with normal models used for G(N ). Motivated by techniques based on first principles for models in discrete time, we propose two models that extend the well known Ricker and Hassel models that satisfy all required properties. The two procedures, the first (G(N ) ⇔ f (r)) based on physical principles and the other on mechanistic principles stemming from the models in discrete time, were interpreted together and provide concepts that allow a better understanding of models for population dynamics. Such concepts, when applied to the musk ox data, have provided interesting insights about the individual interactions starting from population information. A limitation of this study is that only models of a single species were considered. We hope to extend these results to models of various species.
A
Getting G(N ) from f (r)
Using the idea first introduced in [12], the update of the population after a time interval Δt is N (t + Δt) = N (t) + Δt
N
Ii ≡
f (rij )
(39)
j=i
is the interaction strength feeling by the agent i as consequence of its neighbors. Consider dN = ρ(r)dD r as the number of individuals inside a hyper-volume element dD r with position vector r. The density of individuals at r is given by ρ(r). The strength of the interaction felt by the agent i is Ii = dN f (|r − ri |), that is
dD rρ(r)f (|r − ri |).
Ii =
(40)
VD
Then one can show that the per capita growth rate G(N ) = (dN/dt)/N can be computed from the microscopic model by N
1 dD rρ(r)f (|r − ri |), (41) G(N ) = k + N i=1 VD
allowing one to address the macroscopic behavior (per capita growth rate) from the microscopic model f (r).
Getting f (r) from G(N )
B
Let us introduce in this appendix a systematic way to get the microscopic form of the interaction between any two individuals of the population, that is f (r), from the macroscopic properties of the population, given by the per capita growth rate G(N ). To do this, let us introduce the function f˜(r) ≡ rDF −1 f (r) which have an anti-derivative function F˜ (r), that is f˜(r) = dF˜ (r)/dr. One has from the fundamental theorem of calculus that
Rmax
f˜(r)dr = F˜ (Rmax ) − F˜ (Rmin ).
(42)
Rmin
Ri .
(37) The derivative of (42) with respect to N gives
i=1
N Taking Δt → 0, one has dN/dt = i=1 Ri . Using Eq. (7), one gets the ordinary differential equation (ODE) N 1 dN 1 =k+ Ii , N dt N i=1
where
(38)
D1 F D N d F F˜ = dN Ω D ρ0 1 dG(N ) = ρ0 ΩD dN D N =(ΩD /DF )ρ0 r
11
F
(43) (44)
Using the chain rule and the auxiliary variable u ≡ (DF N/(ΩD ρ0 ))1/DF , one can write dF˜ (u) d ˜ du u1−DF f (u) = F (u) = f˜(u) = . dN dr dN Ω D ρ0 Ω D ρ0 (45) Coming back to Eq. (44) and taking u → r, one gets dG(N ) , (46) f (r) = dN D N =μr
where δk,1 is the Kronecker delta, which means that b descendants will be produced in the sites containing exactly one individual, we have the Ricker model: f (Nt ) = bne−Nt /n
Nt /n = bNt e−Nt /n . 1!
To describe spatial aggregation, the negative binomial distribution is often used. It can be described as
F
ΩD D F ρ0 .
with μ ≡ From the above equation one can get the microscopic kind of interaction between the individuals from the macroscopic aspect of the population.
pk =
(Nt /n)k λ Γ(k + λ) Γ(λ) Γ(k + 1) (λ + Nt /n)k+λ
(50)
where Γ(x) is the gamma function and λ is a positive parameter related to the degree of aggregation of the individuals. If λ is small, there will be empty sites and consequently individuals are C The site based framework more more concentrated in clusters. For λ → ∞ we have In this appendix we summarize the results pre- Poisson distribution and if λ = 1 we obtain the sented in [8]. In this article, a systematic method geometric distribution given by k −(k+1) to deduce the Ricker and Hassel models (among Nt Nt . 1+ pk = others) by first principles is discussed. n n Consider n discrete sites where N individuals live. We define pk as the proportion of sites con- If φ(k) = bδk,1 , we have, using equation (50): taining k individuals. The reproduction depends bλ1+λ Nt only on the amount of individuals in the site. φ(k) (51) Nt+1 = λ+1 , is the interaction function, defined as the expected λ + Nnt number of individuals emerging from the site with k individuals. The expected population at time t + 1 which is the Hassel model. Another interaction function is (Nt+1 ) given that there are Nt individuals at time t is given by (52) φ(k) = bkC k−1 ∞ pk φ(k) (47) with 0 ≤ C ≤ 1. The interaction function (49) Nt+1 = f (Nt ) ≡ n is a particular case of (52) when C → 0. When k=0 C → 1, we have the emergence of bk individuals To get the Ricker model, it is sufficient to asat each site of k individuals in the previous gensume that individuals are uniformly distributed. eration, the maximum number possible. We can Making n → ∞ keeping density Nt /n constant, interpret 1 − C ≡ s as related to the struggle bethe number of individuals in a given site follows tween individuals. Combining Eqs. (48) and (52) a Poisson distribution with mean Nt /n, so that we have k (−Nt /n) /k! From (47): pk = (Nt /n) e (53) Nt+1 = bNt e−sNt /n . ∞ (Nt /n)k Similarly, using Eqs. (50) and (52): φ(k). (48) Nt+1 = f (Nt ) = ne−Nt /n k! k=0 bλ1+λ Nt Nt+1 = (54) λ+1 Given φ(k), substituting in Eq. (48) we obtain a λ + s Nnt model for the dynamics in discrete time for Nt . For These are the models of Ricker and Hassel with the interaction function given by a specific parameter associated with the degree of (49) competition, respectively. φ(k) = bδk,1 12
References [1] Allee WC (1931) Animal Aggregations: A Study in General Sociology (University of Chicago Press, Chicago).
[12] Mombach JC, Lemke N, Bodmann BE, Idiart MAP (2002) A mean-field theory of cellular growth. EPL (Europhysics Letters) 59:923.
[13] d’Onofrio A (2009) Fractal growth of tumors and other cellular populations: Linking the [2] Courchamp F, Berec L, Gascoigne J (2008) mechanistic to the phenomenological modelAllee effects in ecology and conservation (Oxing and vice versa. Chaos, Solitons & Fractals ford University Press, Oxford, New York). 41:875–880. [3] Smith R, et al. (2014) Programmed allee effect [14] Ribeiro, Fabiano L. Nascimento K (2014) A in bacteria causes a tradeoff between populaone dimensional model of population growth. tion spread and survival. Proceedings of the Submitted to Physica A. National Academy of Sciences 111:1969–1974. [15] Ribeiro FL (2015) A non-phenomenological [4] Vercken E, Kramer A, Tobin P, Drake J (2011) model of competition and cooperation to exCritical patch size generated by allee effect in plain population growth behaviors. Bulletin of gypsy moth, lymantria dispar (l.). Ecology letmathematical biology 77:409–433. ters 14:179–186. [16] Cross SS (1997) Fractals in pathology. The [5] Brashares JS, Werner JR, Sinclair A (2010) Journal of pathology 182:1–8. Social ‘meltdown’in the demise of an island endemic: Allee effects and the vancouver island [17] Savage VM, Herman AB, West GB, Leu K (2013) Using fractal geometry and universal marmot. Journal of animal ecology 79:965– growth curves as diagnostics for comparing 973. tumor vasculature and metabolic rate with [6] Boukal DS, Berec L (2002) Single-species modhealthy tissue and for predicting responses to els of the allee effect: extinction boundaries, drug therapies. Discrete and continuous dysex ratios and mate encounters. Journal of namical systems. Series B 18. Theoretical Biology 218:375–394. [18] Guiot C, Degiorgis PG, Delsanto PP, Gabriele [7] Ricker WE (1954) Stock and recruitment. P, Deisboeck TS (2003) Does tumor growth Journal of the Fisheries Board of Canada follow a “universal law”? Journal of theoreti11:559–623. cal biology 225:147–151. [8] Br¨ annstr¨om ˚ A, Sumpter DJ (2005) The role of [19] Baish JW, Jain RK (2000) Fractals and cancompetition and clustering in population dycer. Cancer research 60:3683–3688. namics. Proceedings of the Royal Society B: [20] Tsyganov M, et al. (1999) The mechanism Biological Sciences 272:2065–2072. of fractal-like structure formation by bacte[9] Hassell M (1975) Density-dependence in rial populations. Journal of biological physics single-species populations. The Journal of an25:165–176. imal ecology pp 283–295. [21] Nonnenmacher TF, Losa GA, Weibel ER [10] Beverton R, Holt S (1993) On the Dynam(2013) Fractals in biology and medicine ics of Exploited Fish Populations, Chapman & (Birkh¨ auser). Hall fish and fisheries series (Springer Nether[22] Klein K, Maier T, Hirschfeld-Warneken VC, lands). Spatz JP (2013) Marker-free phenotyping of [11] Sanchez A, Gore J (2013) Feedback between tumor cells by fractal analysis of reflection inpopulation and evolutionary dynamics deterterference contrast microscopy images. Nano mines the fate of social microbial populations. letters 13:5474–5479. 13
[23] Falconer K (2004) Fractal Geometry: Mathematical Foundations and Applications (Wiley).
[35] Tener JS (1965) Muskoxen in Canada: a biological and taxonomic review (Department of Northern Affairs and National Resources, Canadian Wildlife Service).
[24] Cabella BCT, Martinez AS, Ribeiro F (2011) Data collapse, scaling functions, and analytical solutions of generalized growth models. Physical Review E 83:061902.
[36] Hubert BA (1977) Estimated productivity of muskox on truelove lowland. Truelove Lowland, Devon Island, Canada: A High Arctic Ecosystem pp 467–491.
[25] Ribeiro FL (2014) A non-phenomenological model to explain population growth behaviors. arXiv preprint arXiv:1402.3676.
[37] Pearce CM (1991) Mapping muskox habitat in the canadian high arctic with spot satellite data. Arctic pp 49–57.
[26] Ribeiro F, Cabella BCT, Martinez AS (2014) Richards-like two species population dynamics model. Theory in Biosciences 133:135–143. [27] Tsoularis A, Wallace J logistic growth models. sciences 179:21–55.
(2002) Analysis of Mathematical bio-
[28] Turchin P (2003) Complex population dynamics: a theoretical/empirical synthesis (Princeton University Press) Vol. 35. [29] Aguirre P, Gonz´alez-Olivares E, S´aez E (2009) Three limit cycles in a leslie-gower predatorprey model with additive allee effect. SIAM Journal on Applied Mathematics 69:1244– 1262. [30] Stephens PA, Sutherland WJ (1999) Consequences of the allee effect for behaviour, ecology and conservation. Trends in Ecology & Evolution 14:401–405. [31] Berec L, Angulo E, Courchamp F (2007) Multiple allee effects and population management. Trends in Ecology & Evolution 22:185–191. [32] Aguirre P, Gonz´alez-Olivares E, S´aez E (2009) Two limit cycles in a leslie-gower predator– prey model with additive allee effect. Nonlinear Analysis: Real World Applications 10:1401–1416. [33] Gregory SD, Bradshaw CJ, Brook BW, Courchamp F (2010) Limited evidence for the demographic allee effect from numerous species across taxa. Ecology 91:2151–2161. [34] Liddle AR (2007) Information criteria for astrophysical model selection. Monthly Notices of the Royal Astronomical Society: Letters 377:L74–L78. 14