A two-sex branching process with oscillations: application to predator-prey systems
aa r X i v : . [ q - b i o . P E ] J a n A two-sex branching process with oscillations: applicationto predator-prey systems
Cristina Guti´errez ∗† Carmen Minuesa ∗‡ January 29, 2021
Abstract
A two-type two-sex branching process is introduced with the aim of describing theinteraction of predator and prey populations with sexual reproduction and promiscuousmating. In each generation and in each species the total number of individuals whichmate and produce offspring is controlled by a binomial distribution with size given bythis number of individuals and probability of success depending on the density of preysper predator. The resulting model enables us to depict the typical cyclic behaviour ofpredator-prey systems under some mild assumptions on the shape of the function thatcharacterises the probability of survival of the previous binomial distribution. We presentsome basic results about fixation and extinction of both species as well as conditions forthe coexistence of both of them. We also analyse the suitability of the process to modelreal ecosystems comparing our model with a real dataset.
Keywords: predator-prey model; two-sex branching process; oscillations; promiscuous mating;extinction; coexistence; density dependence.
MSC:
Recently, the first stochastic process to model the interplay of predator and prey populationswith sexual reproduction was introduced in Guti´errez and Minuesa (2020). Precisely, the modelis a two-type two-sex controlled branching process with a promiscuous mating system and wherethe control of the populations depends on the current number of individuals of each type inthe ecosystem. The model has a recurrent dynamics in each generation where we distinguishthree phases: reproduction, control and mating. First, during the reproduction stage couples ∗ Both authors contributed equally to this work. † Department of Mathematics, University of Extremadura, 10071, C´aceres, Spain. E-mail address: [email protected] . ORCID: 0000-0003-1348-748X. ‡ Department of Mathematics, Autonomous University of Madrid, 28049, Madrid, Spain. E-mail address: [email protected] . ORCID: 0000-0002-8858-3145.
1f each species produce offspring. Next, during the control stage, predators may die due tothe lack of preys and preys can be killed by predators. We also allow the possibility of deathof individuals by other causes (hunting, existence of other predators, other natural causes,etc.) during this phase. The last stage is the mating phase where females and males of eachspecies mate before giving birth to their children during the reproduction phase at the followinggeneration. Particularly, we consider a promiscuous mating system, that is, where a male maymate with several females, but whenever there is some male in the population each femalemates with only one of them. Examples of this kind of mating system in the nature can befound in Norris and Schilt (1988), F¨urtbauer et al. (2011), Balme et al. (2013), Lifjeld et al.(2019), Wightman et al. (2019) or Lee et al. (2019). In this setting, we establish conditions forthe extinction or coexistence of both species and the fixation of one of them. The reader is alsoreferred to Guti´errez and Minuesa (2020) for further motivation and background informationon the literature on predator-prey models.Although this paper constitutes the first approach to the problem of modelling predator-preysystems with sexual reproduction, the model presents some drawbacks. Namely, the probabilityof survival of each individual in both populations depends on the sizes of both species in absoluteterms. Given that the prey-to-predator ratio is known to be a useful predictor of predationrate and prey growth rate (see Vucetich et al. (2011)), a more realistic situation is to considerthat the survival of each individual depends on this density or proportion of preys per predatorand its relation with some parameter µ . Intuitively, this parameter represents the proportionof preys per predator that enables the populations to stay stable and it plays an important rolein the description of the behaviour of the process regarding the coexistence of both species, aswe state below.While this modification retains the generation-by-generation dynamics described above, itleads to a new model that is the focus of this work. An additional advantage of this new processversus the model in Guti´errez and Minuesa (2020) is that it enables us to model the cyclicbehaviour observed in many predator-prey systems. We highlight that this process extends thedynamics of other well-known predator-prey models defined via ordinary differential equations(ODEs) and to that end, we examine the similarities between these models and our predator-prey process. We also illustrate these facts through several simulated examples and show thesuitability of the model to mimic the evolution of real predator-prey systems. Motivated bythe impact of predator-prey systems on human life (see, for instance, Karachle et al. (2016),Boulˆetreau et al. (2018), or Ohlberger et al. (2019)) and the possible applicability of our resultsfor the preservation of species, we study the behaviour of our model in relation to the extinctionproblem. More precisely, we analyse the possible fixation of each species and we providesufficient conditions for the coexistence of both the predator and prey populations.Apart from this introduction, the paper is organised in 6 sections and one appendix. InSection 2 we provide the formal definition of the model and an intuitive interpretation of theassumptions. In Section 3 we study the fixation of each species and the possibility of theultimate extinction of the whole predator-prey system. Section 4 is devoted to the analysis ofthe coexistence of both species. In Section 5 we illustrate the versatility of our process to modelreal world predator-prey systems. We summarise the main results of this work in Section 6.The proofs of all the results are collected in a final appendix to ease the readability of the2aper.In the following, all the random variables (r.v.s) are defined on the same probability space(Ω , A , P ). Moreover, we write N = N ∪ { } , and let A denote the indicator function of theset A . In this section, we are concerned with a two-type two-sex and density-dependent branchingprocess and its application to predator-prey systems. This model evolves as a three-stagemodel of reproduction, control and mating in each generation. We also consider a promiscuousmating as the mating system. We shall start with the formal definition and interpretation ofthe model.Let us denote Z n and e Z n the number of predator couples and prey couples at the n -thgeneration, respectively, and assume that the population starts with an arbitrary number ofpredator and prey couples: ( Z , e Z ) = ( z, ˜ z ) ∈ N . Let us fix a generation n ∈ N . The dynamics of the three phases in the next generationsand the interpretation of the variables involved in the model is as follows. First, in the repro-duction phase , every couple of each species produces offspring independent of the others andin accordance with the same probability law. Let us denote ξ ni the number of offspring of the i -th predator couple at generation n while ˜ ξ ni denotes the number of offspring of the i -th preycouple at generation n , and let us assume that the common probability distributions may varybetween the different species, but it remains constant over the generations for each species.Therefore, { ξ ni : n ∈ N , i ∈ N } and { ˜ ξ ni : n ∈ N , i ∈ N } are two independent families ofr.v.s and the r.v.s. within each family are independent and identically distributed non-negativeand integer valued. At the end of this phase, the sum of all the offspring of each species givesus the total number of predators, T n +1 , and the total number of preys, e T n +1 , in the followinggeneration n + 1: ( T n +1 , e T n +1 ) = Z n X i =1 ξ ni , e Z n X i =1 ˜ ξ ni , (1)where the empty sums are defined as 0.The reproduction stage is followed by the control phase , where the number of predatorsand preys that are available for reproduction might not coincide with the total number ofindividuals due to several reasons (hunting, lack of food supply, their capture by predators,reproduction disability, etc.) and therefore they might be reduced. Thus, if there are T n +1 = t predators and e T n +1 = ˜ t preys in the population, the r.v.s φ n +1 ( t, ˜ t ) and ˜ φ n +1 ( t, ˜ t ) denote thenumber of predators and preys that survive and are able to reproduce in generation n + 1,respectively. Then, we assume that { φ n ( t, ˜ t ) : n, t, ˜ t ∈ N } and { ˜ φ n ( t, ˜ t ) : n, t, ˜ t ∈ N } are twoindependent sequences of independent and non-negative integer valued r.v.s. We also considerthat the survival of each predator (prey) is independent of the generation and of the survival3f the remaining predators (preys), and that the probability of survival is the same for all theindividuals in the same species. Consequently, it is natural to assume binomial distributions forthe control variables. Precisely, if there are t ∈ N predators and ˜ t ∈ N preys in the populationat generation n + 1, then we assume that the distribution of the variable φ n +1 ( t, ˜ t ) is binomialwith size t and probability of success s (˜ t/t ), whereas the one of ˜ φ n +1 ( t, ˜ t ) is also binomial withsize ˜ t and probability ˜ s (˜ t/t ), with s, ˜ s : [0 , ∞ ) → [0 ,
1] continuous functions, that is, for t, ˜ t ∈ N , φ n +1 ( t, ˜ t ) ∼ B ( t, s (˜ t/t )) , ˜ φ n +1 ( t, ˜ t ) ∼ B (˜ t, ˜ s (˜ t/t )) . We note that s (˜ t/t ) and ˜ s (˜ t/t ) represent the probabilities that a predator and a prey survive, re-spectively, given that there are t predators and ˜ t preys in the population, and these probabilitiesdepend on the density of preys per predator, ˜ t/t .At the end of this control phase, there are F n +1 females and M n +1 males within the survivorpredator population, and e F n +1 females and f M n +1 males within the survivor prey populationat generation n + 1 and all of them are able to reproduce. We write α and ˜ α , with 0 <α, ˜ α <
1, to denote the probability that a survivor predator and prey is female, respectively.Thus, if there are k predators and ˜ k preys that survive the control phase, then the randomvector ( F n +1 , M n +1 ) follows a multinomial distribution with parameters k and ( α, − α ), and( e F n +1 , f M n +1 ) follows a multinomial distribution of parameters ˜ k and ( ˜ α, − ˜ α ). Formally,conditional on { φ n +1 ( T n +1 , e T n +1 ) = k, ˜ φ n +1 ( T n +1 , e T n +1 ) = ˜ k } ,( F n +1 , M n +1 ) ∼ M ( k ; ( α, − α )) , ( e F n +1 , f M n +1 ) ∼ M (˜ k ; ( ˜ α, − ˜ α )) . The last step is the mating phase , where female and male predators and preys mate in orderto form couples at generation n + 1, Z n +1 and e Z n +1 . The total number of couples of each speciesis determined by means of a promiscuous mating function ( f ( x, y ) = x min { , y } ) dependingon the number of females and males of each species at the current generation:( Z n +1 , e Z n +1 ) = (cid:16) F n +1 min { , M n +1 } , e F n +1 min { , f M n +1 } (cid:17) . The process { ( Z n , e Z n ) } n ∈ N defined above is named predator-prey two-sex and density-dependent branching process (PP-2SDDBP). Note that by the definition of the model it isnot difficult to verify that the number of predator and prey couples in a certain generation onlydepends on the total number of predator and prey couples in the previous generation. Thus, thebivariate sequence { ( Z n , e Z n ) } n ∈ N is a discrete time homogeneous Markov chain whose statesare two-dimensional vectors with non-negative integer coordinates. Moreover, since the emptysum is assumed to be zero, if there is no couple of one of the species in some generation, fromthis generation on, couples and individuals of that species no longer exist. That is, if Z n = 0( e Z n = 0, resp.) for some n ∈ N , then T k = 0 ( ˜ T k = 0, resp.) and Z k = 0 ( e Z k = 0, resp.) for all k ≥ n . Bearing this in mind, it is not difficult to deduce that (0 ,
0) is an absorbing state andevery non-null state is transient.Finally, let us introduce the notation for the distribution of the offspring variables. Wewrite p k = P ( ξ = k ), and ˜ p k = P ( ˜ ξ = k ), for k ∈ N . The distributions p = { p k } k ∈ N p = { ˜ p k } k ∈ N are called offspring distribution or reproduction law of the predator and preypopulations , respectively. In order to avoid trivialities, we assume that p + p + p < p + ˜ p + ˜ p <
1. Moreover, we assume that these distributions have finite and positive meansand variances, which are denoted m and σ , respectively, for the predators, and ˜ m and ˜ σ forthe preys.A natural question that arises in this setting is what conditions the functions s ( · ) and ˜ s ( · )should satisfy so that the model can describe realistic situations observed in nature and whichare also covered in other predator-prey models, such as, for instance:(S1) The limited appetite of predators (see Arditi et al. (1978)).(S2) The survival of predators in absence of preys due to alternative food supplies or theirartificial feeding (see Adamu (2018) or Prasad et al. (2013)).(S3) The fact that preys have a chance to hide away from predators (see Gal et al. (2015)).(S4) The obvious fact that animals die for natural reasons, that is, not all predators die becauseof the lack of preys and not all preys are killed by predators.(S5) The typical cyclic behaviour in this type of predator-prey systems given by the fact thatwhen the proportion of preys per predator is small (big) the number of preys and predatorsshould drop (arise) in the next generations.To this end, henceforth we assume that these functions satisfy the next conditions for certainconstants 0 < ρ < ρ <
1, 0 < ˜ ρ < ˜ ρ < µ > s ( · ) and ˜ s ( · ) are strictly increasing continuous functions.(C2) lim x →∞ s ( x ) = ρ and lim x →∞ ˜ s ( x ) = ˜ ρ .(C3) lim x → s ( x ) = ρ and lim x → ˜ s ( x ) = ˜ ρ .(C4) The r.v. φ n ( k,
0) follows a binomial distribution with parameters k and ρ , and ˜ φ n ( k,
0) =0 a.s., for each n ∈ N , and k ∈ N .(C5) The r.v. ˜ φ n (0 , ˜ k ) follows a binomial distribution with parameters ˜ k and ˜ ρ , and φ n (0 , ˜ k ) =0 a.s., for each n ∈ N , and ˜ k ∈ N .(C6) s ( µ ) = 1 αm and ˜ s ( µ ) = 1˜ α e m .The PP-2SDDBP { ( Z n , e Z n ) } n ∈ N satisfying conditions (C1)-(C6) is called predator-prey two-sex and density-dependent branching process with oscillations (PP-2SDDBPO). The intuitiveinterpretation of those conditions is now explained.The monotonicity property (condition (C1)) of s ( · ) and ˜ s ( · ) is a natural assumption sinceif the ratio of the preys to the predators is very large, then there are enough resources for thepredators to survive and many preys also survive whenever we assume the limited appetite5f predators (situation (S1)) together with the possibility that preys escape from predators(situation (S3)). Consequently, the probability of survival of each predator and prey should behigh. The greater this ratio becomes, the greater the probability of survival of each predatorand prey. This occurs because the available resources of each predator -in relative terms-increases so does the probability of predator survival. Moreover, each prey has a less chance tobe captured by predators, and as a result the probability that a prey survives arises.Condition (C2) ((C3), resp.) means that when the number of preys is much greater (smaller)than the number of predators in some generation, the probability of survival of predatorsand preys converge to certain constants ρ and ˜ ρ ( ρ and ˜ ρ ). Intuitively, ρ is the survivalprobability of predators when there is no prey and ˜ ρ is the survival probability of preys whenthere is no predator. Moreover: • The condition ρ > • The fact that preys could escape from predators, situation (S3), is described by ˜ ρ > • The assumption ρ < • Similarly, the constraint ˜ ρ < µ interpreted as the necessaryprey density per predator so that both populations stay stable. Suppose that there are t predators and ˜ t preys in the population in certain generation. When the proportion of preys perpredator is less than prey density needed to maintain the predator population, that is, ˜ t/t < µ ,then there is a lack of preys to feed the predators and some of them may die of starvation,thus, the expected number of predators at the next generation should be less than the currentpopulation size. The same behaviour should occur for the prey population because the chanceto be captured by the predators is greater as a result of the large number of predators. Then, inthis case, the number of individuals of both populations should decrease. On the other hand,when the proportion of preys per predator is greater than prey density needed to maintainthe predator population around the same population size, that is, ˜ t/t > µ , then the predatorshave enough preys to feed themselves and survive and consequently, the expected number ofpredators or preys at the next generation should be greater than the current number of thoseindividuals in the ecosystem.The above facts are described through condition (C6). For a better understanding of themeaning this condition we establish the following proposition. The proof is analogous to theone of Proposition 7 (i) in Guti´errez and Minuesa (2020) and it is therefore omitted. Let us6rst introduce the following notation: G n = σ ( T l , ˜ T l : l = 1 , . . . , n ) , n ∈ N . Proposition 2.1.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO and let us fix n ∈ N . Then:(i) The conditional expectations of the number of individuals of each species are E [ T n +1 |G n ] = αmT n s ( e T n /T n ) h − (cid:0) − s ( e T n /T n )(1 − α ) (cid:1) T n − i { T n > } ,E [ e T n +1 |G n ] = ˜ α e m e T n ˜ s ( e T n /T n ) h − (cid:0) − ˜ s ( e T n /T n )(1 − ˜ α ) (cid:1) e T n − i { T n > } + ˜ α e m e T n ˜ ρ h − (cid:0) − ˜ ρ (1 − ˜ α ) (cid:1) e T n − i { T n =0 } . (ii) The conditional variances of the number of individuals of each species satisfy: V ar [ T n +1 |G n ] ≤ h α m T n s ( e T n /T n ) (1 − s ( e T n /T n )(1 − α )) T n − + ( σ + m ) αT n s ( e T n /T n ) i { T n > } ,V ar [ e T n +1 |G n ] ≤ α e m e T n h ˜ s ( e T n /T n ) (1 − ˜ s ( e T n /T n ) + ˜ α ˜ s ( e T n /T n )) e T n − { T n > } + ˜ ρ (1 − ˜ ρ + ˜ α ˜ ρ ) e T n − { T n =0 } i + (˜ σ + e m ) ˜ α e T n h ˜ s ( e T n /T n ) { T n > } + ˜ ρ { T n =0 } i . For an intuitive explanation of condition (C6) let us analyse Proposition 2.1 (i) . Let usfocus on the predator population and consider that both populations have enough number ofindividuals at generation n to avoid extinction at the next generation. In particular, T n > e T n >
1. Since ρ >
0, then 1 − s ( e T n /T n )(1 − α ) ≤ − ρ (1 − α ) <
1. Consequently, if T n becomes very large, then the term T n (cid:0) − s ( e T n /T n )(1 − α ) (cid:1) T n − is negligible and essentially, E [ T n +1 |G n ] ≈ αmT n s ( e T n /T n ) . Now, if the proportion of preys per predator e T n /T n is equal to µ , then it is natural to expectthat the number of predators in the generation n + 1 is roughly the same that in the currentgeneration n , and this can be achieved by setting s ( µ ) = 1 / ( αm ). Under this condition, if e T n /T n > µ , then the expected number of individuals is greater than the current number ofpredators in the population, whereas if e T n /T n < µ , then there is a drop in the expected numberof individuals. The same arguments provide the interpretation for the prey population.Notice that (C1)-(C3) together with (C6) imply: ρ αm < < ρ αm and ˜ ρ ˜ α e m < < ˜ ρ ˜ α e m, (2)7nd therefore αm > α e m >
1. For this reason, throughout the paper we assume that theparameters of the model satisfy (2).
Remark 2.1.
We provide two examples of functions s ( · ) and ˜ s ( · ) satisfying (C1)-(C3). Forinstance, for the function s ( · ) , we could take f ( x ) = ( ρ − ρ )(1 − a − x ) + ρ , with a > , or f ( x ) = ρ x b + ρ x b + 1 , with b > . In order to those functions satisfy condition (C6), we should take a = (cid:18) ( ρ − ρ ) αmρ αm − (cid:19) /µ , and b = log((1 − ρ αm ) / ( ρ αm − µ ) , whenever b > . Note that a > always holds under condition (2) . In this subsection, we present two simulated examples to illustrate different behaviours thatcan be observed in the model. We consider a PP-2SDDBPO with Poisson distributions forboth the predator and prey reproduction laws. We simulated the first n generations of thismodel starting with Z = 10 predator couples and e Z = 10 prey couples. We take µ = 2 inboth examples. The remaining values of the parameters are shown in Table 1. n α m ρ ρ ˜ α e m ˜ ρ ˜ ρ Example 1 75 0.42 3 0.45 0.9 0.49 4 0.25 0.6Example 2 3000 0.49 4 0.25 0.6 0.42 3 0.45 0.9Table 1: Values of the parameters of the models in the simulated examples.In Example 1, we have ˜ ρ ˜ α e m = 1 . > ρ αm = 1 . µ before its fast increase.A theoretical justification for these results is given in Theorem 4.2 in Section 4.8
20 40 60 + + + Generations N u m be r o f c oup l e s PredatorsPreys (in decens) 0 20 40 60 + + + + + Generations T o t a l o f i nd i v i dua l s be f o r e t he c on t r o l pha s e PredatorsPreys (in decens) 0 20 40 60
Generations R a t i o o f p r e ys t o p r eda t o r s be f o r e c on t r o l Figure 1:
Example 1 . Left: evolution of the number of predator couples (solid line) and preycouples (dashed-dotted line). Centre: evolution of the total number of predators (solid line)and of the total number of preys (dashed-dotted line) before the control phase. Right: evolutionof the ratio of the total of preys to the total of predators before the control phase (black line).Horizontal dashed-dotted line represents the value of µ . Generations N u m be r o f c oup l e s PredatorsPreys 0 500 1000 1500 2000 2500 3000
Generations T o t a l o f i nd i v i dua l s be f o r e t he c on t r o l pha s e PredatorsPreys 0 500 1000 1500 2000 2500 3000
Generations R a t i o o f p r e ys t o p r eda t o r s be f o r e c on t r o l Figure 2:
Example 2.
Left: evolution of the number of predator couples (solid line) andprey couples (dashed-dotted line). Centre: evolution of the total number of predators (solidline) and of the total number of preys (dashed-dotted line) before the control phase. Right:evolution of the ratio of the total of preys to the total of predators before the control phase(black line). Horizontal dashed-dotted line represents the value of µ .In Example 2, we show the case when ρ αm = 1 . > ˜ ρ ˜ α e m = 1 . µ beforethey hit the value 0 at generation 2785. Despite the existence of the mating phase, we now make a comparison between our model andthe equivalent deterministic model defined by the Lotka-Volterra (predator-prey) equations. Tothat end, we analyse the evolution of the expected number of individuals of each species thatsurvive in presence of the other one. We shall start with the prey population.Let us assume that there are x predators and ˜ x preys in the ecosystem, then the meannumber of preys that are alive after the control phase is ˜ x ˜ s (˜ x/x ) = ˜ x − ˜ x ˜ d (˜ x/x ), with ˜ s ( · ) =1 − ˜ d ( · ), and ˜ d ( · ) denoting the probability that a prey does not survive. Thus, the term˜ x corresponds to the growth of the number of preys in the Lotka-Volterra equations whereas˜ x ˜ d (˜ x/x ) -expected number of preys died during the control- is the counterpart of the interactionterm in the Lotka-Volterra equations. We also note that growth rate for the preys in theLotka-Volterra equations is introduced in our model via the reproduction phase of this species.Moreover, if at some generation there is no predator, then by Proposition 2.1 E [ e T n +1 | T n = 0 , e T n = ˜ x ] = ˜ α e m ˜ x ˜ ρ h − (cid:0) − ˜ ρ (1 − ˜ α ) (cid:1) ˜ x − i , ˜ x ∈ N , which implies an exponential growth of the prey population as occurs in predator-prey systemof ODEs.Next, we examine the mean number of predators which are alive after the control phase.This quantity is xs (˜ x/x ), which corresponds to the interaction term in the Lotka-Volterraequations. The term corresponding to the death of predators by other causes different from thelack of preys is introduced in our model by the consideration that ρ < predation rate (PR) and the kill rate (KR). The former represents the proportion of the preypopulation killed by the predators and the latter is the number of preys killed per predatorduring a unit of time. Since the present model is a discrete-time stochastic process, thesequantities refer to expected quantities in each generation. Thus, the equivalent of the predatorrate in our model is the probability of survival ˜ d (˜ x/x ), whereas the counterpart of the kill rateis ˜ d (˜ x/x )˜ x/x , that is, the expected number of preys caught per predator during the control10hase. We emphasise the fact that these quantities satisfies the common relationship P R = KR · P dP r , where
P d stands for the number of predators and
P r for the number of preys.Finally, we remark that on the one hand, the fact that (0 ,
0) is an absorbing state is theequivalent in our model of the fixed point (0 ,
0) for the Lotka-Volterra system, which representsthe ultimate extinction of both populations. On the other hand, the parameter µ in thePP-2SDDBPO (recall that µ is the prey-to-predator ratio that ensures that both populationsstay stable) is the counterpart of the second steady state in the Lotka-Volterra equations,where both populations support their current size (see Hirsch and Smale (1974), pp. 258-263for further details). Moreover, the introduction of the parameter µ together with the shapeof the functions s ( · ) and ˜ s ( · ) enables us to model the fluctuations usually observed in thesedynamical systems, as we showed in Example 2 in the previous subsection. In this section we examine the evolution of the species in relation to their fixation and extinction.For this purpose, we first consider the following sets: { Z n → , e Z n → } , termed extinctionof predators and prey populations; { Z n → ∞ , e Z n → ∞} , coexistence of both species; { Z n →∞ , e Z n → } , predator population fixation; and { Z n → , e Z n → ∞} , prey population fixation.In a first result we establish the usual property of branching processes known as theextinction-explosion. The proof of this result is similar to the proof of Proposition 8 inGuti´errez and Minuesa (2020) and it is therefore omitted. Proposition 3.1.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO. Then:(i) P (lim inf n →∞ ( Z n , e Z n ) = ( k, ˜ k )) = 0 , and P (lim sup n →∞ ( Z n , e Z n ) = ( k, ˜ k )) = 0 , for each ( k, ˜ k ) ∈ N \{ (0 , } .(ii) P ( Z n → , e Z n →
0) + P ( Z n → ∞ , e Z n → ∞ ) + P ( Z n → ∞ , e Z n →
0) + P ( Z n → , e Z n →∞ ) = 1 . Taking into account this result, we examine the behaviour of both species in relation to theirextinction. To that end, given i, j >
0, henceforth we write P ( i,j ) ( · ) = P ( ·| ( Z , e Z ) = ( i, j )).First, we study the probability of fixation of each species. In the fixation events, from somegeneration on, the survivor species behaves as the two-sex branching process with randomcontrol on the total number of individuals introduced in Section 2 in Guti´errez and Minuesa(2020). The control variables of the corresponding processes follow binomial distributions withconstant probability of success ρ , in the case of the predator fixation, and with probability˜ ρ in the case of the prey fixation. Thus, bearing in mind condition (2) the following resultimmediately holds applying Theorems 1 and 2 (ii) (b) in the aforementioned paper.11 roposition 3.2. Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO. For any initial values i, j > :(i) P ( i,j ) ( Z n → ∞ , e Z n →
0) = 0 .(ii) P ( i,j ) ( Z n → , e Z n → ∞ ) > . This result establishes that the fixation of the predator population is not possible becausethe mean growth rate of this species in absence of the prey population is too small ( ρ αm < ρ ˜ α ˜ m > Remark 3.1.
We emphasize that Proposition 3.2 (i) is not in contradiction with situation (S2)because (S2) refers to the survival of each predator individual, while Proposition 3.2 (i) concernsthe predator population. That is, it is possible that one predator survives in absence of preys,but in the long run, if there is no prey in the ecosystem the growth rate of predators is notenough to prevent the population from the extinction.
The following result on the ultimate extinction of the predator-prey system follows imme-diately by Proposition 3.2.
Corollary 3.1.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO. Then, P ( i,j ) ( Z n → , e Z n → < forany initial values i, j > . We recall that the extinction of both populations is always possible due to many reasonssuch as the possibility that no individual of both populations survives during the control phase,or that all the survivors of both species are of the same sex, which makes impossible to formnew couples. In addition, there is a positive probability that the predator and prey couplesproduce no offspring if p > p > We now examine the possibility of the coexistence of both species. Although until now, wehave focussed on the fate of the total number of couples, Z n and e Z n , of each species, the resultson coexistence in terms of the total number of individuals of each species before the controlphase, T n and e T n , are equivalent in view of Lemma 1 in appendix. Thus, for our conveniencein the proofs of the results in this section we establish the theorems in terms of T n and e T n .Our results show that the coexistence of both species strongly depends on the parameter µ .For that reason, we writeΩ = n lim sup n →∞ e T n T n ≤ µ o ∪ n lim inf n →∞ e T n T n > µ o ∪ n lim inf n →∞ e T n T n ≤ µ < lim sup n →∞ e T n T n o . (3)12n the first theorem of this section we establish that the survival of none species is possiblein the first event of (3). As a consequence, the coexistence is also impossible for any initialnumber of individuals and regardless the values of the parameters. The idea is that if thenumber of preys per predator is not greater than µ from some generation on, then, on the onehand, the predator population will die out owing to the lack of resources to survive, and, onthe other hand, more preys will be captured given the needs of the large number of predatorscompared to the prey population. Theorem 4.1.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO. For any initial values i, j > : P ( i,j ) n lim sup n →∞ e T n T n ≤ µ o ∩ { T n → ∞} ! = 0 ,P ( i,j ) n lim sup n →∞ e T n T n ≤ µ o ∩ { e T n → ∞} ! = 0 . In particular, P ( i,j ) n lim sup n →∞ e T n T n ≤ µ o ∩ { T n → ∞ , e T n → ∞} ! = 0 . In our second result, we provide sufficient conditions for the coexistence of both species. Thisstates that when ˜ ρ ˜ α e m (the maximum mean growth rate of female preys) is greater than ρ αm (the maximum mean growth rate of female predators), then the coexistence of both species ispossible in the second event in (3). This result was illustrated in Example 1 in Subsection 2.1. Theorem 4.2.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO. If any of the following two conditionshold: (i) ρ αm < ˜ ρ ˜ α e m , (ii) ρ αm = ˜ ρ ˜ α e m , and lim inf x →∞ (˜ s ( x ) e m ˜ α ) / ( s ( x ) mα ) > ,then for any initial values i, j > : P ( i,j ) (cid:18)n lim inf n →∞ e T n T n > µ o ∩ { T n → ∞ , e T n → ∞} (cid:19) > . In particular, P ( i,j ) ( { T n → ∞ , e T n → ∞} ) > . The previous results do not cover the possibility of coexistence in the case ρ αm > ˜ ρ ˜ α e m .Our numerical studies indicate that if the prey fixation does no occur at early generations ofthe process, then the process enters in an oscillation phase until the number of preys is zero,and then the ultimate extinction of the entire system occurs. In particular, in that phase the13atio e T n /T n fluctuates around µ (see Figure 2 in Example 2). We therefore expect that thecoexistence is impossible on both events { lim inf n →∞ e T n /T n > µ } and { lim inf n →∞ e T n /T n ≤ µ < lim sup n →∞ e T n /T n } if ρ αm > ˜ ρ ˜ α e m . Moreover, our simulation analysis suggests that theoscillation stage cannot occur forever, regardless the value of the parameters. We formalizethose ideas in the following conjecture. Conjecture 4.1.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBPO. For any initial values i, j > : P ( i,j ) (cid:18)n lim inf n →∞ e T n T n ≤ µ < lim sup n →∞ e T n T n o ∩ { T n → ∞ , e T n → ∞} (cid:19) = 0 . Moreover, if ρ αm > ˜ ρ ˜ α e m , then for any initial values i, j > : P ( i,j ) (cid:18)n lim inf n →∞ e T n T n > µ o ∩ { T n → ∞ , e T n → ∞} (cid:19) = 0 , and in particular, P ( i,j ) (cid:0) T n → ∞ , e T n → ∞ (cid:1) = 0 . To support this conjecture we performed a simulation of 10 processes until generation 5 · ,with the same parameters as in Example 2 in Subsection 2.1. We estimated the probabilityof coexistence of both species by using the proportion of processes, among the 10 simulatedones, where both species have survived. We show its evolution over the generations in Figure 3,where we observe how that probability decreases to zero as the number of generations increases. . . . . Genenerations P r obab ili t y o f c oe x i s t en c e Figure 3: Evolution of the probability P ( T n > , e T n >
0) based on the simulation of 10 processes following a simulated PP-2SDDBPO with ρ αm = 1 . > ˜ ρ ˜ α e m = 1 . µ = 2. 14 Modelling of the population of wolves and moose inIsle Royale
In this section, we evaluate the suitability of the proposed model to describe the evolutionof the number of wolves and moose studied in Isle Royale National Park over 61 years (seeVucetich and Peterson (2012)). The Isle Royal National Park is an international biospherereserve located in Lake Superior (Michigan, USA). The wolf and moose populations in thisecosystem are well known since their study started back in 1959. Indeed, it is one of thepredator-prey system with the longest continuous study in the world. In the following, we showthat this real dataset can be easily mimicked using the model presented in this paper. To thatend, we generated a path of a PP-2SDDBPO where the values of the parameters are reasonablebearing in mind what one can observe in populations of wolves and moose. We remark thatour aim here is to illustrate the adjustment and usefulness of our process to model the realworld situations, but making inference on the parameters of the model exceeds the scope ofthis paper.First, we note that the reproduction, control and mating phases described in the Section 2faithfully adjust to the evolution of the populations of wolves and moose. This is owing tothe fact that both species have a sexual reproduction with promiscuous mating where not allfemales are able to reproduce (see, for example, Ballenberghe and Miquelle (1993) for the moosepopulation or Peterson et al. (2002) for the wolf population).Next, we describe the choice of the parameter values and their relationship with the popu-lations. We shall start with the reproduction laws. On the one hand, given that several authorsagreed that there is an unbalanced sex ratio in wolves, with more males than females in thepopulation (see Cowan (1947), Mech (1975) or Stenlund (1955)), we considered α = 0 .
35 assex ration for the predator. On the other hand, we took ˜ α = 0 .
49 as the sex ratio for the preysmotivated by the fact some studies claimed that there is no statistically significant differencebetween the population of males and of females in moose populations (see Peek et al. (1976)or Shubin and Yazan (1959)).Regarding the reproduction, since a typical wolf litter consists of four to seven pups, for themean of the predator reproduction law we fixed m = 4. Similarly, we took e m = 2 . ρ = ˜ ρ = 0 .
05, both of them very close to 0. The first choice isdue to the fact that the Isle Royale is an isolated place where the moose is the main food sourcefor the wolves. Thus, their survival is very difficult in absence of moose (see Peterson and Page(1988)). The latter is because the probability that the moose could avoid predators is small ifthe proportion of wolves per prey goes to infinity. For the probabilities ρ and ˜ ρ , we considered ρ = 0 . ρ = 0 .
93. The choice ρ = 0 . ρ = 0 .
93 isjustified by several facts. First, because wolves are their main (or even their unique) predator;second, because moose hunting is prohibited in the island, and third, because the possibilityof immigration and emigration is negligible (see Vucetich and Peterson (2004)). We note thatthe choice of the parameters is reasonable in the context of Isle Royale biosphere, but they also15atisfy conditions (2) and ρ αm < ˜ ρ ˜ α e m .Finally, given that our simulation studies show that the proportion of preys per predatoroscillates around the parameter µ in the present model, we therefore considered that µ = 35 is aplausible value as illustrated in Figure 4 (right). For the PP-2SDDBPO with these parameterswe simulated 61 generations starting with the same number of individuals as the real dataset,that is, 20 predators and 538 preys. We remark that neither the initial number of couples northeir offspring are known for this real predator-prey system. Moreover, we fixed the function f ( · ) in Remark 2.1 to define the probabilities of survival ( s ( · ) and ˜ s ( · )).First, in Figure 4 (left) we plot the total number of wolves and the total number of dozens ofmoose observed in the Isle Royale National Park from 1959 to 2019. Next, in Figure 4 (centre)we show the evolution of the total number of predators and the total number of dozens of preys ofthe simulated path in each generation. To obtain this last path, we simulated processes followingthe model described above, and we found that it is easy to get paths with a similar behaviour tothe one showed in Figure 4. Note that the total number of individuals which is observed in thepopulation is previous to the control phase, that is, they correspond to the values of the variables { ( T n , e T n ) : n = 1 , . . . , } in our model. Some of these individuals will participate in themating phase but, initially, we do not have any knowledge about how many of them. The totalnumber of the individuals after the control, that is, { ( φ n ( T n , e T n ) , ˜ φ n ( T n , e T n )) : n = 1 , . . . , } ,is therefore unknown. Generations N u m be r o f i nd i v i dua l s Moose/12Wolves 0 10 20 30 40 50 60
Generations N u m be r o f i nd i v i dua l s Preys/12Predators 0 10 20 30 40 50 60
Generations P r opo r t i on o f p r e y pe r p r eda t o r Real d ata Simulation
Figure 4: Left: evolution of the number of wolves (dashed line) and moose (solid line) from1959 to 2019 in the real dataset. Centre: evolution of the number of predators (dashed line)and of the number of preys (solid line) previous to the control phase over 61 generations in thesimulated process. Right: proportion of preys per predator over the generations for both thereal data (solid line) and the simulated path (dashed line). The horizontal line represents thevalue µ = 35.By comparing both graphs one can observe a similar behaviour. Overall, the total numberof dozens of preys is always greater than the total number of predators over the generations forboth the real and the simulated dataset. Moreover, the ratio between the number of individualsin both populations remains oscillating around some constant value over the generations. Thisfact is an indicator that the model can capture correctly the evolution of the proportion of16redators per prey, as can be seen in Figure 4 (right). In relation to the cyclic behaviour ofour model, one can notice that stages of growth in the number of preys are followed by stagesof growth in the predator population and vice-versa, as it occurs with the real data. Moreover,when the number of preys drops, the same tendency is appreciated in the number of predatorsin the next generations.It is worth to mention that the real dataset presents a dramatic drop in moose populationin the mid 1990s, from a total of 2398 in 1996 to 900 in 1997. Some authors claimed that thissharp drop in the real data might have been caused by the severe climatic conditions during thatperiod (see Vucetich and Peterson (2004) or Montgomery et al. (2013)). We remark that ourmodel is not able to capture this type of extreme changes because the offspring distributions donot vary over the generations, and then, the drops and rises of our model are more moderated.Nevertheless, one could think of including a varying environment for the reproduction laws todescribe these situations, however, we do not analyse this fact in depth because it is beyondthe scope of this paper. The aim of this work is to model the generation-by-generation evolution of populations ofpredators and preys as well as their interaction when both of them have sexual reproduction. Forthat purpose, a two-type two-sex branching process with a control in the number of individualsof each species is introduced in the present paper.As usual in two-sex branching processes theory, the reproduction and mating phases aretaking into account in the definition of the model. Additionally, a control stage is also consid-ered to model the relationship between the species. Those three stages are repeated in eachgeneration. The description for a given generation is the following: first, in the reproductionphase, a random number of offspring is produced by couples of each species. Next, due to theinteraction between predators and preys, some of these individuals could die. This is modelledin the control phase by means of binomial distributions where the probability of success foreach species depends on the density of preys per predator. To determine how many survivorindividuals are females and males multinomial distributions are considered. In this final matingphase, couples are formed via a promiscuous mating in both species. These couples produceoffspring in the next generation and this three-stage cycle is repeated again.We have proved some results concerning with the fate of the species in the population(extinction, fixation and coexistence). In particular, we have established that for any initialnumber of couples of both species the probability of fixation of the predator population isnull due to the fact that the mean growth rate of this species is less than one in absence ofpreys. Additionally, the prey population has a positive probability of fixation because its meangrowth rate is greater than one in absence of predators, and as a consequence, the probabilityof extinction is less than one.The study of the possibility of the coexistence is more complex due to its dependence onthe parameter µ . This parameter represents the necessary prey density per predator so thatboth populations remain stable. Specifically, we have established that the coexistence is not17ossible if the supremum of the density of preys per predator is not above µ in the long run.Furthermore, when the infimum of such a density is eventually greater than µ , we provide asufficient condition for the coexistence to have a positive probability. Based on a simulatedstudy, we have also conjectured that, on the one hand, the oscillation phase cannot last foreveron the coexistence event. On the other hand, when the maximum mean growth rate of femalepredators is greater than the maximum mean growth rate of female preys, and both populationsare still alive, the process fluctuates during long time until its final extinction and then thecoexistence of both species is not possible.The main novelty of this model regarding to the one introduced in Guti´errez and Minuesa(2020) for population with sexual reproduction and also with respect to other predator-preymodels in the context of branching process (see Alsmeyer (1993) or Coffey and B¨uhler (1991)) isthat under certain conditions this process presents the typical fluctuations or oscillations whichappear in nature in the majority of predator-prey systems (see, for example, Gilpin (1973)).In those oscillations we observe periods with a large number of predators and preys which arefollowed by periods of a smaller number of predators and preys. Two elements are the key ofthis model to obtain this behaviour: first, the introduction of the density of preys per predatorin the survival function in the control phase of each species and second, the introduction ofthe parameter µ representing the necessary proportion of preys per predator to maintain bothpopulations balanced.As a final conclusion, we emphasize that this model has shown to be appropriate to fit realpredator-prey ecosystems. The fact of considering sexual reproduction between the individu-als of species results particularly useful to model the interaction of mammals populations ofpredators and preys. Funding
This research was funded by the Ministerio de Ciencia e Innovaci´on [grant PID2019-108211GBI00/AEI/10.13039/501100011033].
Appendix
Lemma 1.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBP. Then: { T n → ∞ , e T n → ∞} = { Z n → ∞ , e Z n → ∞} a.s. Proof of Theorem 4.1
We shall prove the first part. The second one is obtained with the same arguments. We followthe same steps as in the proof of Proposition 10 in Guti´errez and Minuesa (2020); we providethe details for the readers convenience. First, note that ( lim sup n →∞ e T n T n ≤ µ ) ⊆ ∞ [ N =1 ( sup n ≥ N e T n T n ≤ µ ) , N > P ( i,j ) ( sup n ≥ N e T n T n ≤ µ ) \ n T n → ∞ o! = 0 . (4)Next, observe that for all n ∈ N , using Proposition 2.1 (ii) we have E [ T n +1 |G n ] ≤ αmT n s ( e T n /T n ) ≤ T n a.s. on { e T n /T n ≤ µ } . (5)Now, let us fix N >
0, and introduce the sequence of r.v.s { X n } n ∈ N defined as: X n = ( T N + n , if N + n ≤ τ ( µ ) ,T τ ( µ ) , if N + n > τ ( µ ) , n ∈ N , where τ ( µ ) is the stopping time: τ ( µ ) = ( ∞ , if sup n ≥ N e T n T n ≤ µ, min (cid:8) n ≥ N : e T n T n > µ (cid:9) , otherwise . If we prove that { X n } n ∈ N is a non-negative super-martingale with respect to {G N + n } n ∈ N ,then applying the martingale convergence theorem, we have that the sequence { X n } n ∈ N con-verges almost surely to a non-negative and finite limit X ∞ , where X ∞ = ( lim n →∞ T n , if sup n ≥ N e T n T n ≤ µ,T τ ( µ ) , otherwise.As a consequence, (4) is obtained and the proof of the first part finishes.To prove that { X n } n ∈ N is a super-martingale with respect to {G N + n } n ∈ N , we first note thatthe variable X n is G N + n -measurable for any n ∈ N . Next, we show that E [ X n +1 |G N + n ] ≤ X n a.s. for each n ∈ N .Let us fix n ∈ N . Then, we have that if e T N + k T N + k ≤ µ , for every k = 0 , . . . , n , then the stoppingtime satisfies τ ( µ ) ≥ N + n + 1, and using (5) we get E [ X n +1 |G N + n ] = E [ T N + n +1 |G N + n ] ≤ T N + n = X n a.s. on n \ k =0 ( e T N + k T N + k ≤ µ ) . Now, if there exists k ∈ { , . . . , n } satisfying e T N T N ≤ µ, . . . , e T N + k − T N + k − ≤ µ and e T N + k T N + k > µ , then forthe stopping time we have τ ( µ ) = N + k ≤ N + n < N + n + 1, and consequently, for theconditional expectation we get E [ X n +1 |G N + n ] = E [ T τ ( µ ) |G N + n ] = T τ ( µ ) = X n a.s. on k − \ l =0 ( e T N + l T N + l ≤ µ ) \ ( e T N + k T N + k > µ ) . n ∈ N , the stopping time is τ ( µ ) = N < N + n + 1 if e T N T N > µ , and, consequently, E [ X n +1 |G N + n ] = E [ T N |G N + n ] = T τ ( µ ) = X n a.s. on ( e T N T N > µ ) . Combining all the above, if A n = ∩ nk =0 (cid:8) e T N + k T N + k ≤ µ } is a measurable set with respect to the σ -algebra G N + n , then we get E [ X n +1 |G N + n ] = E [ X n +1 |G N + n ] A n + E [ X n +1 |G N + n ] A cn ≤ T N + n A n + T τ ( µ ) A cn = X n a.s. (cid:3) In order to prove Theorem 4.2, we make use of the following easy-to-prove lemma.
Lemma 2.
Let { ( Z n , e Z n ) } n ∈ N be a PP-2SDDBP. If p + p + p < and ˜ p + ˜ p + ˜ p < ,then the sets { ( i,
0) : i > } , { (0 , j ) : j > } and { ( i, j ) : i, j > } are classes of communicatingstates and each state leads to the state (0 , . Furthermore, the process can move from the lastset to the others in one step. Proof of Theorem 4.2
The proof follows similar arguments to those in Theorem 3 of Guti´errez and Minuesa (2020).Let us fix N ∈ N . We shall prove that we can take I , J ∈ N sufficiently large so that P ( i,j ) (cid:18)n inf n ≥ N e T n T n > µ o ∩ { T n → ∞ , e T n → ∞} (cid:19) > , for each i, j ∈ N , i ≥ I , j ≥ J . (6)If this holds, then the result also holds for any i, j ∈ N . Indeed, if i, j ∈ N satisfy that either i < I or j < J , then the existence of n ∈ N such that P ( Z n = I , e Z n = J | Z = i, e Z = j ) > , is guaranteed by Lemma 2, and using the Markov property we have P ( i,j ) (cid:18)n lim inf n →∞ e T n T n > µ o ∩ { T n → ∞ , e T n → ∞} (cid:19) ≥≥ P ( i,j ) (cid:18)n inf n ≥ N + n e T n T n > µ o ∩ { T n → ∞ , e T n → ∞} (cid:19) ≥ P ( I ,J ) (cid:16) inf n ≥ N e T n T n > µ, T n → ∞ , e T n → ∞ (cid:17) · P ( i,j ) ( Z n = I , e Z n = J ) > . Hence, it only remains to prove (6). We shall start with the case ˜ ρ ˜ α e m > ρ αm >
1. Then,we can choose ε , ε > < ρ αm − ε = η < ˜ ρ ˜ α e m − ε = η . p + p + p < p + ˜ p + ˜ p <
1, we can also take r , r ≥ p r > p r >
0, and then, η = r − > η = r − >
1. Moreover, sincelim x →∞ s ( x ) = ρ , and lim x →∞ ˜ s ( x ) = ˜ ρ , there exists M > µ such that | αmρ − αms ( x ) | ≤ ε / | ˜ α e m ˜ ρ − ˜ α e m ˜ s ( x ) | ≤ ε /
4, for every x ≥ M .We now consider the next events: A = { η Z < T , ˜ η e Z < e T } and A n = { η T n < T n +1 , η e T n < e T n +1 , M T n ≤ e T n } , n ∈ N . Then, we have that P ( i,j ) (cid:16) inf n ≥ N e T n T n > µ, T n → ∞ , e T n → ∞ (cid:17) ≥ P ( i,j ) (cid:16) inf n ∈ N e T n T n > µ, T n → ∞ , e T n → ∞ (cid:17) ≥ P ( i,j ) ( ∩ ∞ n =0 A n )= lim n →∞ P ( i,j ) ( ∩ nl =0 A l )= lim n →∞ P ( i,j ) ( A ) n Y l =1 P (cid:16) A l | ∩ l − k =0 A k ∩ { Z = i, e Z = j } (cid:17) . Let us also introduce the sets B l = { ( x, y ) ∈ N : iη η l − < x, j ˜ η η l − < y, M x ≤ y } , for each l ∈ N , and observe that the family ∩ l − k =0 A k ∩ { ( T l , e T l ) = ( x, y ) } ∩ { Z = i, e Z = j } , x, y ∈ N , is a partition of the set ∩ l − k =0 A k ∩ { Z = i, e Z = j } . Then, using the Markov property it is notdifficult to deduce that P (cid:0) A l | ∩ l − k =0 A k ∩ { Z = i, e Z = j } (cid:1) ≥ inf ( x,y ) ∈ B l P (cid:16) A l | ∩ l − k =0 A k ∩ { ( T l , e T l ) = ( x, y ) } ∩ { Z = i, e Z = j } (cid:17) = inf ( x,y ) ∈ B l P (cid:16) A l | ( T l , e T l ) = ( x, y ) (cid:17) = inf ( x,y ) ∈ B l P ( A | ( T , e T ) = ( x, y )) . To finish the proof, we start bounding the probability P ( A c | ( T , e T ) = ( x, y )) from abovewhen ( x, y ) ∈ B l . Observe that if ( x, y ) ∈ B l , then P (cid:0) e T /T < M | ( T , e T ) = ( x, y ) (cid:1) = 0. Next,using Proposition 2.1 (i) we have that for x ∈ N , P (cid:0) T ≤ η T | ( T , e T ) = ( x, y ) (cid:1) == P (cid:16) x (cid:0) ε − αmρ + αms ( y/x ) − αms ( y/x )(1 − s ( y/x )(1 − α )) x − (cid:1) E [ T | ( T , e T ) = ( x, y )] − T (cid:12)(cid:12) ( T , e T ) = ( x, y ) (cid:17) . On the one hand, by the properties of the function s ( · ), we have αms ( y/x )(1 − s ( y/x )(1 − α )) x − ≤ αmρ (1 − ρ (1 − α )) x − , and since 0 < − ρ (1 − α ) <
1, there exists I ∈ N satisfying that αms ( y/x )(1 − s ( y/x )(1 − α )) x − < ε / , for every x ≥ I , y ∈ N . On the other hand, | αmρ − αms ( y/x ) | ≤ ε /
4, for every x, y ∈ N , with y ≥ M x . Thus,combining all the above, for each x ≥ I and y ≥ M x , by Chebyschev’s inequality we have P (cid:0) T ≤ η T | ( T , e T ) = ( x, y ) (cid:1) ≤≤ P (cid:16) x (cid:0) ε − αmρ + αms ( y/x ) − αms ( y/x )(1 − s ( y/x )(1 − α )) x − (cid:1) ≤ | E [ T | ( T , e T ) = ( x, y )] − T | (cid:12)(cid:12) ( T , e T ) = ( x, y ) (cid:17) ≤ P (cid:16) xε ≤ | E [ T | ( T , e T ) = ( x, y )] − T | (cid:12)(cid:12) ( T , e T ) = ( x, y ) (cid:17) ≤ V ar [ T | ( T , e T ) = ( x, y )] ε x . Taking into account Proposition 2.1 (ii) we can choose two constants C > C > x ≥ I and y ≥ M x , P (cid:0) T ≤ η T | ( T , e T ) = ( x, y ) (cid:1) ≤ C x + C (1 − ρ (1 − α )) x . Using the same arguments, we can take J ≥ M I such that for y ≥ J , x ∈ N , with y/x ≥ M , we have that P (cid:0) e T ≤ η e T | ( T , e T ) = ( x, y ) (cid:1) ≤≤ P (cid:16) y (cid:0) ε − ˜ α e m ˜ ρ + ˜ α e m ˜ s ( y/x ) − ˜ α e m ˜ s ( y/x )(1 − ˜ s ( y/x ) + ˜ α ˜ s ( y/x )) y − (cid:1) ≤ | E [ e T | ( T , e T ) = ( x, y )] − e T | (cid:12)(cid:12) ( T , e T ) = ( x, y ) (cid:17) ≤ P (cid:16) yε ≤ | E [ e T | ( T , e T ) = ( x, y )] − T | (cid:12)(cid:12) ( T , e T ) = ( x, y ) (cid:17) ≤ C y + C (1 − (1 − ˜ α ) ˜ ρ ) y , for some constants C > C > P ( i,j ) ( A ) from below observe that P ( i,j ) ( A ) = P ( i,j ) ( T > η Z , e T > ˜ η e Z ) ≥ P ( T = r i | Z = i ) · P ( e T = r j | e Z = j ) = p ir ˜ p jr > , ⌊ η i ⌋ + 1 ≤ r i and ⌊ ˜ η i ⌋ + 1 ≤ r i , for any i ∈ N , and ⌊·⌋ denotes theinteger part function.Finally, if i ≥ I and j ≥ J , then since η , ˜ η , η , η >
1, we have iη η l − ≥ I and j ˜ η η l − ≥ J , for all l ∈ N . Now, it is not difficult to verify that P ( i,j ) ( T n → ∞ , e T n → ∞ ) ≥ p ir ˜ p jr ∞ Y l =1 inf ( x,y ) ∈ B l P ( A | ( T , e T ) = ( x, y )) ≥ p ir ˜ p jr · ∞ Y l =1 inf ( x,y ) ∈ B l (cid:18) − C x − C (1 − (1 − α ) ρ ) x − C y − C (1 − (1 − ˜ α ) ˜ ρ ) y (cid:19) ≥ p ir ˜ p jr · ∞ Y l =1 (cid:18) − C iη η l − − C (1 − (1 − α ) ρ ) iη η l − − C j ˜ η η l − − C (1 − (1 − ˜ α ) ˜ ρ ) j ˜ η η l − (cid:19) > . Now we turn to the case ˜ ρ e m ˜ α = ρ mα , and lim inf x →∞ (˜ s ( x ) e m ˜ α ) / ( s ( x ) mα ) >
1. Then, wecan choose M > µ such that ˜ s ( x ) e m ˜ α > s ( x ) mα for every x ≥ M . Let η = s ( M ) mα > ε = ρ mα − η >
0. Moreover, there exists M ′ > µ such that | ρ mα − s ( x ) mα | ≤ ε / x ≥ M ′ .Next, we take η = ˜ s ( M ) e m ˜ α , and ε = ˜ ρ e m ˜ α − η >
0, which satisfies η > η > M > M ′ such that | ˜ ρ e m ˜ α − ˜ s ( x ) e m ˜ α | ≤ ε /
4, for every x ≥ M . Theproof continues with the same arguments as before and the result follows. (cid:3) References
H. A. Adamu. Mathematical analysis of predator-prey model with two preys and one predator.
International Journal of Engineering and Applied Sciences , 5(11):17–23, 2018.G. Alsmeyer. On the Galton-Watson predator-prey process.
The Annals of Applied Probability ,3(1):198–211, 1993.R. Arditi, J. M. Abillon, and J. Vieira Da Silva. A predator-prey model with satiation andintraspecific competition.
Ecological Modelling , 5(3):173–191, 1978.V. V. Ballenberghe and D. G. Miquelle. Mating in moose: timing, behavior, and male accesspatterns.
Canadian Journal of Zoology , 71:1687–1690, 1993.G. A Balme, A. Batchelor, N. de W. Britz, G. Seymour, M. Grover, L. Hes, D. W. Macdon-ald, and L. T. B. Hunter. Reproductive success of female leopards Panthera pardus: theimportance of top-down processes.
Mammal Review , 43(3):221–237, 2013.23. Boulˆetreau, A. Gaillagot, L. Carry, S. T´etard, E. De Oliveira, and F. Santoul. Adult Atlanticsalmon have a new freshwater predator.
PLOS ONE , 13(4):1–12, 2018.J. Coffey and W. J. B¨uhler. The Galton-Watson Predator-Prey process.
Journal of AppliedProbability , 28(1):9–16, 1991.I. M. Cowan. The timber wolf in the Rocky Mountain National Parks of Canada.
CanadianJournal of Research , 25(5):139–174, 1947.A. A. Danilkin and A. A. Ulitin. A review of moose fertility in Russia.
Alces , 34(2):459–466,1998.I. F¨urtbauer, M. Heistermann, O. Sch¨ulke, and J. Ostner. Concealed fertility and extendedfemale sexuality in a non-human primate (Macaca assamensis).
PLoS ONE , 6(8):e23105,2011.S. Gal, S. Alpern, and J. Casas. Prey should hide more randomly when a predator attacksmore persistently.
Journal of The Royal Society Interface , 12(113):20150861, 2015.M. E. Gilpin. Do Hares Eat Lynx?
The American Naturalist , 107(957):727–730, 1973.C. Guti´errez and C. Minuesa. A predator-prey two-sex branching process.
Mathematics , 8(9):1408, 2020. doi: 10.3390/math8091408.M. W. Hirsch and S. Smale.
Differential Equations, Dynamical Systems, and Linear Algebra ,volume 60 of
Pure and Applied Mathematics . Academic Press, Inc., 1974.P. Karachle, A. Angelidis, G. Apostolopoulos, D. Ayas, M. Ballesteros, Bonnici C., M. M.Brodersen, L. Castriota, N. Chalari, J. M. Cottalorda, F. Crocetta, A. Deidun, Z. Dodo,A. Dogrammatzi, D. Jakov, F. Fiorentino, O. Gonulal, J. G. Harmelin, G. Insacco, andA. Zenetos. New Mediterranean biodiversity records (March 2016).
Mediterranean MarineScience , 17(1), 2016.N. S. Lee, N. L. Goodwin, K. E. Freitas, and A. K. Beery. Affiliation, aggression, and selectivityof peer relationships in meadow and prairie voles.
Frontiers in Behavioral Neuroscience , 13:52, 2019.J. T. Lifjeld, J. Gohli, T. Albrecht, E. Garcia-Del-Rey, L. E. Johannessen, O. Kleven, P. Z.Marki, T. C. Omotoriogun, M. Rowe, and A. Johnsen. Evolution of female promiscuity inPasserides songbirds.
BMC Evolutionary Biology , 19(1):169, 2019.L. D. Mech. Disproportionate sex ratios in wolf pups.
Journal Wildlife Manage , 39:737–740,1975.R. A. Montgomery, J. A. Vucetich, R. O. Peterson, G. J. Roloff, and K. F. Millenbach. Theinfluence of winter severity, predation and senescence on moose habitat use.
Journal ofAnimal Ecology , 82(2):301–309, 2013. 24. S. Norris and C. R. Schilt. Cooperative societies in three-dimensional space: On the originsof aggregations, flocks, and schools, with special reference to dolphins and fish.
Ethology andSociobiology , 9:149–179, 1988.J. Ohlberger, D. E. Schindler, E. J. Ward, T. E. Walsworth, and T. E. Essington. Resurgenceof an apex marine predator and the decline in prey body size.
Proceedings of the NationalAcademy of Sciences , 116(52):26682–26689, 2019. doi: 10.1073/pnas.1910930116.J. M. Peek, D. L. Urich, and R. J. Mackie. Moose habitat selection and relationships to forestmanagement in Northeastern Minnesota.
Wildlife Monographs , (48):3–65, 1976.R. O. Peterson and R. E. Page. The rise and fall of Isle Royale wolves, 1975ˆa € “1986. Journalof Mammalogy , (69):89–99, 1988.R. O. Peterson, A. K. Jacobs, T. D. Drummer, L. D. Mech, and D. W. Smith. Leadership behav-ior in relation to dominance and reproductive status in gray wolves, Canis lupus.
CanadianJournal of Zoology , (80):1405–1412, 2002.B. S. R. V. Prasad, M. Banerjee, and P. D. N. Srinivasu. Dynamics of additional food providedpredatorˆa € “prey system with mutually interfering predators. Mathematical Biosciences , 246(1):176–190, 2013.G. G. Shubin and P. Yazan. Experiment in the organization and conduct of a moose tradefarm.
Proc. Pechora-Ilych State Nature Preserve , 7:288–324, 1959.M. H. Stenlund. A field study of the timber wolf (Canis lupus) on the Superior National Forest,Minnesota.
Minnesota Department of Conservation Technical Bulletin , 4, 1955.J. A. Vucetich and R. O. Peterson. The influence of top–down, bottom-up and abiotic factorson the moose (Alces alces) population of Isle Royale.
Proceedings. Biological sciences , 271:183–189, 2004.J. A. Vucetich and R. O. Peterson. The population biology of Isle Royale wolves and moose:an overview. , 2012.J. A. Vucetich, M. Hebblewhite, D. W. Smith, and R. O. Peterson. Predicting prey popula-tion dynamics from kill rate, predation rate and predator-prey ratios in three wolf-ungulatesystems.
Journal of Animal Ecology , 80:1236–1245, 2011.P. H. Wightman, J. C. Kilgo, M. Vukovich, J. R. Cantrell, C. R. Ruth, B. S. Cohen, M. J.Chamberlain, and B. A. Collier. Gobbling chronology of eastern wild turkeys in SouthCarolina.