A predator-prey two-sex branching process

In this paper, we introduce a two-sex controlled branching model to describe the interaction between predator and prey populations with sexual reproduction. This process is a two-type branching process, where the first type corresponds to the predator population and the second one to the prey population. While each population is described via a two-sex branching model, the interaction and survival of both groups is modelled through control functions depending on the current number of individuals of each type in the ecosystem. We provide necessary and sufficient conditions for the ultimate extinction of both species, the fixation of one of the species and the coexistence of both of them. Moreover, the description of the present predator-prey two-sex branching process on the fixation events can be performed in terms of the behaviour of a one-type two-sex branching process with a random control on the number of individuals, which is also introduced and analysed.


Introduction
Predator-prey models have been widely studied in the literature since the introduction of the first model in [32]. This process aimed at modelling the trophic interactions between both species and as later in [44], non-linear differential equations were used to describe the interaction of a predatorprey dynamic system. Since then, many modifications and new models have been defined trying to adapt the peculiarities observed in the real world as faithfully as possible. Some of the most recent papers on this topic are, for example, [28], [31], [14], [47], and [40], which deal with theoretical questions, or [11] and [20], which focus on real world settings.
All those papers use deterministic models based on ordinary differential equations (ODEs) to model the predator-prey dynamics. However, the interaction between both species can be seen as a stochastic system. In this sense, [23] conducts an interesting study about the dynamics of deterministic and stochastic models for a predator-prey system, where the predator species suffers from a parasitic infection. The deterministic model is an ODE model while the stochastic model is derived by means of continuous time Markov chains. The theory of branching processes is used to estimate the probability of a disease outbreak and the probability of the prey species invasion.
In the context of branching processes several publications have tackled predator-prey system modelling. We stand out the pioneer work of [22] where two predator-prey models in continuous time are considered: a host-parasite system and a predator-prey process for which the predator birth rate is not directly associated with a prey death. A further study of the process in [22] is performed in [39]. In [7], a two-type branching process in discrete time is introduced to describe two populations, predators and preys, living in the same ecosystem. In the model, the evolution of the predator population is independent of the population size of preys, and the number of preys at each generation is given by the number of prey offspring minus the number of preys that have been captured and killed by the predators. The analysis of this process continues in [1] and it focusses on the study of necessary and sufficient conditions for the fixation of both populations. Later, [6] analyses a counterpart of the predator-prey model in [7] in the continuous-time setting by making use of linear birth-death processes. More recently, in [4] a branching model is used to analyse and compare the influence of the habitat lost, poaching or drop of preys in tiger populations. Moreover, it is important to mention the work of [10] where a branching random walk is used to describe predator-prey populations.
However, none of the previous models takes into account the fact that many animal species have a sexual reproduction. This fact plays a key role in the evolution of the species because their development not only depends on the number of individuals of the species, but also on the total number of females and males and the type of mating of the populations. Many real cases have been reported on the interaction between preys and predators where both populations are formed by females and males which mate and procreate by means of sexual reproduction. For example, [34] studied populations of sea lions and penguins as predators and preys, respectively, and in [35], a social network of giraffe populations is studied bearing in mind the presence of lions which preferentially prey on giraffe calves.
Our aim in this paper is to introduce a predator-prey stochastic process by using the two-sex branching processes theory (see [8]) in order to fill this niche in the literature on predator-prey models. Therefore, the novelty of the present paper is to consider a two-type two-sex discrete-time branching process in order to model the predator-prey interaction of populations of females and males with sexual reproduction. Moreover, we focus on situations where each female mates with only one male, whenever there are males in the population, but a same male could mate with more than one female. This type of mating is called promiscuous mating and we can observe it in many examples in the nature (see, for example, [36], [13], [2], [30], [46] or [27]).
The two-sex branching process allows us to describe the generation-by-generation evolution of the number of individuals of a predator-prey system in certain environment. The definition of the model is based on the following assumptions: in each generation, females and males of each species mate and form couples by promiscuous mating, and each of those couples is assumed to give birth to some number of individuals. Nevertheless, the survival of all these individuals to form couples and reproduce is constrained due to the interaction between both species. Thus, some preys could be captured and killed by the predators to feed themselves and some predators could die due to lack of food supplies. As a result, the couples of each species of the following generation will be formed from the females and males that have survived.
The second aim of this work is to study how the number of individuals of each species evolves over successive generations. We examine conditions for one of the species -predator or prey-to become extinct or to have a positive probability of survival. Moreover, we are interested in studying sufficient conditions for the coexistence of both species assuming two different scenarios. First, we analyse the destiny of the species assuming that a predator population with limitless appetite is introduced and regarded as an invasive species in a geographically isolated area where its only food supply is the prey population. Second, we assume that predators with limited appetite could survive without any prey due to the presence of other food resources. This consideration on the number of preys consumed per predator will give the preys a chance to survive although they might coexist with a large predator population. This paper is divided into 7 sections. Apart from this introduction, a controlled two-sex branching process is introduced in Section 2. This model will be useful to facilitate the understanding of the predator-prey two-sex branching model which is presented in Section 3. The basic properties and main results about the conditional moments of the variables involved in the predator-prey two-sex branching model are studied in Section 4. In Sections 5 and 6, we provide conditions for the ultimate extinction of the entire predator-prey system, for the coexistence of both species or for the fixation of one of them; in Section 5 we focus on the case when both species live together in an isolated environment, whereas in Section 6 we consider that they might coexist with other animal species in a non-isolated ecosystem. To conclude, a discussion about the main results obtained is reported in Section 7. In order to facilitate the reading, the proofs of the results are gathered in some appendices.

A controlled two-sex branching process
Before providing the definition of the predator-prey two-sex branching process, we begin by introducing a simpler branching process that will play an important role in the analysis developed in the following sections. The model is a generalization of the two-sex branching process introduced by [8] with the novelty that not all individuals generated by the couples of the previous generation participate in the mating process. Only some selected individuals -males and females-are chosen to formed couples. The selection of those individuals is done through random control functions. Thus, this model represents a combination of the two-sex branching process and the controlled branching process with random control functions on the total number of individuals. First, we provide the probabilistic definition of the model.
Let us consider two independent families of random variables (r.v.s) {ξ ni : n ∈ N 0 , i ∈ N} and {ϕ n (k) : n, k ∈ N 0 }. The former is assumed to be a family of independent and identically distributed (i.i.d.) r.v.s. The latter is assumed to be a family of r.v.s such that {ϕ n (k) : k ∈ N 0 }, n ∈ N 0 , are independent stochastic processes with the same one-dimensional distribution. Let us also consider the sequences {Y n } n∈N 0 and {X n } n∈N defined as: where N 0 = N ∪ {0}, L : N 2 0 → N 0 is a deterministic mating function and (F n+1 ,M n+1 ) is a random vector that follows a multinomial distribution with parameters y and (λ, 1 − λ) conditionally on {ϕ n+1 (X n+1 ) = y}, with 0 < λ < 1. The process {Y n } n∈N 0 is known as two-sex branching process with random control on the total number of individuals (BBPCI). Moreover, the empty sum in (1) is assumed to be 0 and the mating function is assumed to be monotonic and non-decreasing in each argument and it satisfies the following conditions: (A2) L(0, y) = L(x, 0) = 0, for each x, y ∈ N 0 . Intuitively, the variable Y n represents the number of couples at generation n, while the variable X n denotes the total number of individuals -females and males-at generation n. The functions ϕ n (·), n ∈ N 0 , represent a control on the total number of individuals in the population at each generation. Thus, we distinguish three consecutive phases at each generation in this model: control phase, mating phase and reproduction phase. The first phase is a control stage where the number of individuals that participate in the following phase is determined. This number is denoted as ϕ n (X n ) at generation n. Each one of these individuals can be female or male. It will be a female with probability λ and a male with probability 1 − λ. The second phase is the mating phase whereF n females andM n males at generation n mate to produce Y n couples considering the promiscuous mating function as the mating function L, which is defined as L(x, y) = x min{1, y}. The last phase is a reproduction phase where the couples give birth to their offspring. Thus, the variable ξ ni denotes the number of offspring of the i-th couple in the n-th generation. The probability distribution of this variable is named the offspring distribution or reproduction law and it is denoted {q k } k∈N 0 . To avoid trivialities it is assumed to satisfy q 0 + q 1 + q 2 < 1.  We note that although the results in this paper are provided for the promiscuous mating, the majority of them can be easily adapted for these other mating functions.
By the definition of the model, it is not difficult to check that the process {Y n } n∈N 0 is a discrete homogeneous Markov process whose states are non-negative integers. Moreover, if in some generation there are not any couples, that is, Y n = 0 for certain n > 0, and assuming that ϕ 0 (0) = 0 a.s. then, from that generation on, there will be neither individuals nor couples, i.e. X k = 0 and Y k = 0 for all k > n. This implies that the state 0 is absorbing and also the extinction of the population. Similar results about the classification of the states, the extinction and the asymptotic behaviour to those for the controlled branching process (see [16,Chapters 3 and 4]) can be obtained. In the next results we only provide those that will be useful for the proofs of the results in the following sections. The proof of the first proposition is easily obtained with standard procedures and it is omitted. (ii) The classic duality extinction-explosion holds, that is, P (Y n → 0) + P (Y n → ∞) = 1.
An immediate consequence of the previous theorem is that a controlled two-sex branching process with control functions ϕ 0 (k) following binomial distributions with parameters k and 0 < ρ < 1 becomes extinct a.s. if λµρ ≤ 1. This remark will be useful in the proof of the results in Section 6.
(i) For each n ∈ N 0 , the conditional expectation is: where h k (·) denotes the probability generating function (p.g.f.) of ϕ 0 (k).
(ii) If the control function ϕ 0 (k) follows a binomial distribution with parameters k and ρ, then: (a) If f (·) denotes the p.g.f of the variable ξ 01 , then there exist constants 3 Definition of a predator-prey two-sex branching model Having described the controlled two-sex branching process, in this section we introduce a predatorprey two-sex branching process in order to model a predator-prey system. We aim at introducing a model for a biological system where two animal species live together in the same environment, where one of them is the prey and the other one is its natural predator. We focus on the case that both species have sexual reproduction and propose a controlled two-sex branching process to model the evolution of each of them, where the control mechanism is introduced in order to describe their natural interaction. We shall start with the formal definition of the model.
The process defined above enables us to model a predator-prey system with non-overlapping generations. This process evolves as a three-stage procedure of reproduction, control and mating in each generation, where the previous variables have the next interpretation. The variables T n andT n represent the total number of predators and preys, respectively, at generation n, whereas F n and M n (F n andM n ) are the total numbers of progenitor predator females and predator males (prey females and prey males), respectively, at generation n. Moreover, Z n andZ n denote the number of predator couples and prey couples at the n-th generation, respectively. The dynamics of the three phases is described below.
In the reproduction phase, couples of each species produce offspring independent of each others and in accordance with an offspring law. The offspring law may vary for the different species, but it remains constant over the generations for each species. Formally, the number of offspring produced by a couple of each species are represented by sequences of i.i.d. N 0 -valued r.v.s {t ni : n ∈ N 0 , i ∈ N}, and {t ni : n ∈ N 0 , i ∈ N}, where t ni denotes the number of offspring of the i-th predator couple at generation n whilet ni denotes the number of offspring of the i-th prey couple at generation n. The common probability distributions of these variables, p = {p k } k∈N 0 and p = {p k } k∈N 0 , respectively, are called offspring distribution or reproduction law of the predator and prey population, respectively and we assume that they have finite and positive mean and variance, which are denoted m and σ 2 , respectively, for the predators andm andσ 2 for the preys. We recall that we consider that the mating functions L andL are the promiscuous mating function and hence, in order to avoid trivialities, we also assume that p 0 + p 1 + p 2 < 1 andp 0 +p 1 +p 2 < 1. At the end of the reproduction phase at generation n + 1, the sum of all the offspring of each species gives us the total number of predators, T n+1 , and the total number of preys,T n+1 , at this generation.
The reproduction stage is followed by the control phase, where the number of predators and preys could be reduced due to several reasons such as their death because of the hunting, the lack of food supply or their capture by predators. Thus, if there are T n+1 = t predators andT n+1 =t preys in the population, the number of predators and preys which survive are given by the r.v.s φ n+1 (t,t) andφ n+1 (t,t), respectively. Considering that the survival of each predator (prey) is independent of the survival of the remaining predators (preys), and the probability of survival is the same for all the individuals in the same population, then it is natural to assume that the distributions of the variables φ n+1 (t,t) andφ n+1 (t,t) are binomial distributions.
More specifically, φ n+1 (t,t) is assumed to follow a binomial distribution with parameters t and s(t), where the s(t) represents the survival probability of a predator given that there aret preys in the population. The condition on the monotonicity of the function s : R → [0, 1] means that the smaller the number of preys is, the smaller the probability of survival of the predators. Regarding the conditions in (2), we note that ρ 1 > 0 means that the predators could survive although there is no prey in the population because they could find another food source (other prey species, for instance). However, ρ 1 = 0 implies the extinction of the predator population when there are not any preys in the population at some generation. Moreover, ρ 2 < 1 means that, although there are enough preys in the population, the predators could die by several reasons (for example, by hunting or their own predators), whereas ρ 2 = 1 indicates the survival of all the predators when the number of preys goes to infinity.
Analogously, we assume thatφ n+1 (t,t) follows a binomial distribution with parameterst and s(t), withs(t) representing the survival probability of a prey given that there are t predators in the population. The condition on the monotonicity of the functions : R → [0, 1] for the probability of survival of a prey indicates that the greater the number of predators is, the smaller the probability of survival of the preys becomes. The assumptions in (3) also have an intuitive interpretation. The conditionρ 1 > 0 means that the preys could survive although the number of predators goes to infinity because predators have limited appetite. The opposite situation, whenρ 1 = 0, implies the extinction of the prey population when there are a huge number of predators in the population at some generation. In addition,ρ 2 < 1 means that, despite the absence of predators in the ecosystem, the preys could die by another reasons (hunting or another predators). This case is excluded ifρ 2 = 1, which leads to the survival of all preys when there are not any predators in the environment in some generation.
At the end of this control phase, there are F n+1 females and M n+1 males within the survivor predator population at generation n + 1. Thus, if α denotes the probability that a survivor predator is female, then the random vector (F n+1 , M n+1 ) follows a multinomial distribution of parameters y and (α, 1 − α), given that T n+1 = t,T n+1 =t and φ n+1 (t,t) = y. Similarly, there areF n+1 females andM n+1 males within the survivor prey population at generation n + 1, and consequently, (F n+1 ,M n+1 ) follows a multinomial distribution of parametersỹ, and (α, 1 −α), given that T n+1 = t,T n+1 =t andφ n+1 (t,t) =ỹ, and whereα is the probability that a survivor prey is female.
The last step is the mating phase, where the predator and prey couples at generation n+1, Z n+1 andZ n+1 , are determined by means of promiscuous mating functions depending on the number of females and males of each species at the current generation.
Remark 3.1. One can propose several functions s(·) ands(·) satisfying the previous assumptions. For example, for the survival of the predator population and similarly, for the survival of the prey populatioñ

Basic properties of the model
In this section, we establish some basic properties of the process regarding the classification of the states and its main moments. First of all, note that from the definition of the model it is not difficult to deduce that the number of predators and prey couples in a certain generation only depend on the total number of predators and prey couples in the previous generation. Thus, the bivariate sequence {(Z n ,Z n )} n∈N 0 is a discrete time homogeneous Markov chain whose states are two-dimensional vectors with non-negative integer coordinates. Moreover, it is immediate that (0, 0) is an absorbing state taking into account (iii) and (iv) and condition (A2).
In the following easy-to-prove proposition we state some properties of the states associated with the process. (iii) If p 0 + p 1 + p 2 < 1,p 0 +p 1 +p 2 < 1 and 0 = ρ 1 < ρ 2 = 1, then the sets {(0, j) : j > 0} and {(i, j) : i, j > 0} are classes of communicating states and each state leads to the state (0, 0). Furthermore, the states belonging to the second set may move to the other one and to the set {(i, 0) : i > 0} in one step. Finally, the process moves from the last set to the state (0, 0) in one step.
Next, we provide some results concerning the conditional moments of the variables involved in the definition of the process which will be useful in Section 6. Note that from the definition of the model, it is immediate to get the mean and variance of the control variables. Indeed, for n, x,x ∈ N 0 , the conditional expectations of the control variables are and the conditional variances of the control variables are For the next results, we introduce the following notation concerning the σ-algebras generated by the variables involved in the definition of the model. In particular, we denote Then, we have that F n−1 ⊆ G n ⊆ H n , for n ∈ N. Now, the conditional moments of the number of predator and prey individuals can be easily obtained from the definition of the model and hence, the proof of the next proposition is omitted.
Next, we establish some results concerning to the conditional moments of the total number of female and male predators and female and male preys. (ii) The conditional variances of the number of predator and prey females and males given the number of individuals and couples are respectively.
(ii) The conditional variance of the number of predator and prey females and males given the number of couples are respectively.
We note that from this result, it is easy to obtain an upper bound for the conditional expectation of the number of predator and prey couples as stated in the following proposition. Remark 4.1. Since ρ 1 ≤ s(t) ≤ ρ 2 , andρ 1 ≤s(t) ≤ρ 2 , for all t ≥ 0, we obtain upper and lower bounds for the previous conditional expectations, Similar arguments let us get the bounds for the remaining conditional expectations. Now, we provide the conditional moments of the number of predator and prey individuals given the total number of individuals and couples in previous generations. (ii) V ar[T n+1 |G n ] ≤ (σ 2 + m 2 )αT n s(T n ) + 2α 2 m 2 T 2 n s(T n ) 2 (1 − s(T n ) + αs(T n )) Tn−1 .
To conclude this section we establish the following result where we derive the usual property of branching processes known as the extinction-explosion dichotomy.
The sets {Z n → 0,Z n → 0}, {Z n → ∞,Z n → ∞}, {Z n → ∞,Z n → 0} and {Z n → 0,Z n → ∞} are termed extinction of both populations, survival (or coexistence) of both populations, predator population fixation, and prey population fixation, respectively. Moreover, if we denote the extinction and survival of the predator population as {Z n → 0}, and {Z n → ∞}, and the extinction and survival of the prey population as {Z n → 0}, and {Z n → ∞}, then in view of Proposition 4.7 it is immediate that, Remark 4.2. Note that, in the cases of fixation of the predator and the prey, the process behaves as a BBPCI defined (1) in Section 2.

Predator-prey isolated system
In this section we consider an isolated predator-prey system, that is, we assume that both species live together in an isolated area where the prey population constitutes the only food resource for the predators. More specifically, we focus on the case that there is an autochthonous species living in an isolated region and a new (invasive) species is introduced in that ecosystem and it preys on the autochthonous one. The question of how those populations evolve together is tackled in this section and it is of great interest for the preservation of the species in these environments. Several examples of those geographically isolated regions have been reported such as, for instance, Azores Islands (see [33]), Eastern Island (see [41]), Macquarie Island (see [12]), isolated regions in Finland and Northwest Russia (see [43]) or seafloor plateau (see [38]).
In terms of our model, this situation can be expressed as ρ 1 = 0, since the fact that the prey is the only food supply for the predators implies that the probability of survival of any predator is zero if there is no prey in the population. Analogously, we considerρ 1 = 0 which means an unlimited appetite of the predators that implies a null probability of survival for all preys when the number of predators in the population goes to infinity. We also allow ρ 2 ≤ 1 andρ 2 ≤ 1 because both the predators and preys could die due to natural causes although there is no prey or predator, respectively, in the ecosystem.
Notice that, under these assumptions if the number of prey couples at some generation is equal to zero then the prey population becomes extinct forever, that is, ifZ n = 0 for some n > 0, theñ T k = 0 andZ k = 0 for all k > n. The extinction of the prey population also bounds the predator population to disappear due to the fact that s(0) = 0, and consequently, φ k (t, 0) = 0, for all k > n, from which one deduces Z k = 0, and T k+1 = 0, for all k > n. On the other hand, if at some generation n there are no predator couples, then the predator population becomes extinct forever, i.e., T k = 0, and Z k = 0, for all k > n. As a result, and sinces(0) = ρ 2 ,φ k (t, 0) follows a binomial distribution with parameterst and ρ 2 , for all k > n, which means that the prey population behaves as a BBPCI (see Section 2).
Bearing in mind these considerations, we study the fate of both species in the population. To that end, given i, j > 0, we write P (i,j) (·) = P (·|(Z 0 ,Z 0 ) = (i, j)), in the remaining of this paper.

The certain extinction of the predator population
The following result is very natural from the definition of the model and it means that in this kind of populations we cannot have the extinction of the prey population and the survival of the predator population (predator fixation). Recall that the prey population is the only food supply of the predator population, and hence, the extinction of preys dooms the predator population to the extinction.
The following result establishes that in this kind of systems the coexistence of preys and predators is not possible. Intuitively, this is deduced as follows. When the number of predators is too large, the probability of survival of the preys is too low and consequently, we expect a huge drop in the number of preys and as a result, in the number of predators in the next generation. This cycle can be repeated until a generation where the probability of survival of each prey is so negligible that the entire population of preys becomes extinct, so does the predator population in the following generation.
Note that by (5), the previous results and Proposition 4.7 (ii) we obtain the certain extinction of the predator population, that is, P (i,j) (Z n → 0) = 1, for any initial values i, j > 0.

The fixation of the prey population
From Propositions 5.1 and 5.2, we deduce that the predator population becomes extinct almost surely in this model. The question that arises is whether the prey population has a chance to survive depending on its reproductive capacity once the predator population has become extinct. Let us start with the caseρ 2 = 1. On the prey fixation set, the prey population behaves as a standard two-sex branching process (without any kind of control) from one generation on once the predator population has become extinct. Thus, the theory developed in [8] can be applied in this setting and we immediately obtain the following result.
On the other hand, ifρ 2 < 1, on the prey fixation set the prey population behaves as the BBPCI introduced in Section 2 from one generation on. The control variables in this model follow binomial distributions with size equal to the number of preys at the corresponding generations and probabilityρ 2 . In this case, the result is a direct consequence of Theorems 2.1 and 2.2 (ii) (b) provided in the mentioned section.

Predator-prey non-isolated system
In this section we consider the case that 0 < ρ 1 < ρ 2 ≤ 1 and 0 <ρ 1 <ρ 2 < 1. Thus, contrary to Section 5, the predators have a positive probability of survival in absence of the prey population (ρ 1 > 0) due to the availability of other food resources and also a limited appetite (ρ 1 > 0) which allows the prey population to have a positive probability of survival even when the predator population size goes to infinity. Moreover, individuals of the prey population might not reproduce because of the presence of predators, but also for other reasons such as their hunting, diseases, or migratory movements (ρ 2 < 1).

The fixation of the predator and prey populations
In this subsection we study necessary and sufficient conditions for the fixation of each species, that is, for one of the two species (predator or prey) to survive and the other one to become extinct. In the fixation events, the surviving species behaves as the BBPCI introduced in Section 2 from some generation on. The corresponding offspring distribution is the reproduction law of the survivor species and the control functions follows binomial distributions with constant probability of success γ, where γ = ρ 1 in the case of the predator fixation and γ =ρ 2 in the case of the prey fixation. Thus, by using Theorems 2.1 and 2.2 (ii) (b) we have the following result. Proposition 6.1. Let {(Z n ,Z n )} n∈N 0 be a predator-prey two-sex branching process. For any initial values i, j > 0: (i) P (i,j) (Z n → ∞,Z n → 0) > 0 if and only if ρ 1 αm > 1.
Intuitively, this result states that a necessary and sufficient condition for the predator population to have a positive probability of fixation is that the mean number of female predators which survive after the control is greater than one. Alternatively, the second part of the result also states that a necessary and sufficient condition for the prey population to have a positive probability of fixation is that the mean number of female preys which survive after the control is greater than one. Now, taking into account (5) and (6), an immediate consequence of Proposition 6.1 is the following corollary: Corollary 6.1. Let {(Z n ,Z n )} n∈N 0 be a predator-prey two-sex branching process. For any initial values i, j > 0:

The extinction of the population
Intuitively, it is clear that if the mean number of female predators (preys) which survives after the control is less than one then the reproductive capacity of the species is not enough to keep it alive by its own. Therefore, we can establish the following result: Proposition 6.2. Let {(Z n ,Z n )} n∈N 0 be a predator-prey two-sex branching process. For any initial values i, j > 0: Taking into account (5) and (6) and from Proposition 6.2, we deduce the following result on the coexistence of the species: Corollary 6.2. Let {(Z n ,Z n )} n∈N 0 be a predator-prey two-sex branching process. For any initial values i, j > 0, if min{ρ 2 αm,ρ 2αm } ≤ 1, then P (i,j) (Z n → ∞,Z n → ∞) = 0.
We note that there is always a positive probability for the complete extinction of the predatorprey system. This could happens for several reasons: either because all individuals of both species might die during the control phase or because all the survivors of both populations might be of the same sex, which makes impossible to form new couples. If we also allow p 0 > 0 andp 0 > 0, then there is a positive probability that the predator and prey couples produce no offspring.
In the next result, we determine a necessary and sufficient condition for both species to become extinct with probability one, which means the extinction of the entire predator-prey system. Proposition 6.3. Let {(Z n ,Z n )} n∈N 0 be a predator-prey two-sex branching process. For any initial values i, j > 0, P (i,j) (Z n → 0,Z n → 0) = 1 if and only if max{ρ 1 mα,ρ 2mα } ≤ 1.

The predator and prey coexistence
In the following result, we study the possibility of having the coexistence of the predator and prey populations.
Theorem 6.1. Let {(Z n ,Z n )} n∈N 0 be a predator-prey two-sex branching process. For any initial values i, j > 0: In the case of coexistence, both populations grow geometrically at certain rate determined by the parameters of the model (see Figure 1). Note that this kind of growth is the classical behaviour in branching processes including the Galton-Watson predator-prey process (see [1]), but it is not typical in predator-prey systems modelled through ODEs, where periodic cycles are observed. However, there are examples of populations with this exponential growth behaviour (see [43], [42], [37], [45] or [3]). Moreover, our model can be also applied at initial stages of other populations when there is a small number of individuals, such as, for example, in populations of endangered species where the exponential growth is shown (see [26], [21] or [29]).

Discussion
In this paper a controlled two-sex branching process is introduced with the aim of modelling predator-prey interactions in different situations. The generation-by-generation evolution of each species is also studied. We have focussed on the case where both populations have sexual reproduction and females and males form couples via a promiscuous sexual mating system. The process evolves in three stages which are repeated at each generation. First, couples of each species produce a random number of offspring (reproduction phase). Next, some of these individuals die because of the interaction between preys and predators and the remaining ones survive and are able to mate (control phase). This control mechanism is modelled through binomial distributions where for each species the probability of success depends on the number of individuals of the other one. Similarly, once that the number of survivors of each species is known, multinomial distributions are applied to determine the number of females and males among them, where the probability vector is the pair given by the probability that one individual is female and the probability that it is male. The final step at any generation is the mating phase when the couples that are able to produce offspring in the next generation are formed.
We highlight that the main novelty of our model with respect to previous ones is the sexual reproduction among the individuals within each species. This important feature had not been considered in the literature yet, but it is quite common in the nature. We also remark that although we have considered a promiscuous mating, other mating systems could be applied and our results could be adapted to those cases modelling different situations.
To analyse the temporal evolution of the process, we have considered two possible scenarios. In the first one we have assumed that the preys live in an isolated region and a new species, the predator, is introduced. In the second one, both species live together in certain area where predators have different food sources (not only the prey) and a limit appetite; this enables the survival of preys even when there is a large number of predators. On the one hand, our results show that the prey population has two possible behaviours (extinction or explosion) in both settings and this behaviour depends on the mean number of prey females that survives during the control phase. On the other hand, the fate of the predator population depends on the situation. In the first scenario, the predator group becomes extinct almost surely because the prey is its only food supply and the predator has limitless appetite. In the second one, the predator could become extinct or growth indefinitely depending on the mean number of female predators that survive during the natural control of the population.
As mentioned above, the model presented in this paper captures the fact that an increment (drop) in the number of preys implies an increment (drop) in the predator population size, but it does not describe the fluctuation behaviour where a large number of preys and predators is followed after some time by a smaller number of preys and predators (see the classical example of the snowshoe hare and the Canadian lynx in [15]). These type of oscillations are usually caused by fluctuations in the environmental conditions or periodic changes in the reproductive capacity of the species. For instance, in the case of the snowshoe hare, the litter sizes vary between years and females only give birth during the breeding season, which is stimulated by new plants, and it begins around mid-March and runs to August (see [9] or [25]). Then, an appropriate way to introduce this periodic reproductive behaviour in the model is to let the predator offspring distribution and the prey offspring distribution change over the time.
Moreover, even in those populations where an exponential growth is observed, sooner or later the growth stops due to environmental conditions or other causes (for instance, the growth of a population of turtles truncated due to an oil spill is reported in [24]). If we wish to reflect the saturation of the environment -and cut the initial exponential growth-it is necessary to modify our assumptions. For example, we could introduce a carrying capacity parameter in the probability functions s(·) ands(·) or even let them depend on two values simultaneously: the number of the preys and the number of the predators. The introduction of all the aforementioned modifications would lead to a more complex model which is beyond of the scope of this paper and it is left for future research.

Proof of Theorem 2.1
The proof follows the arguments of Theorem 1 in [18]. Since and hence, lim inf n→∞ Y n is finite a.s., from which we get P (Y n → ∞|Y 0 = y) = 0.
Proof of Theorem 2.2 where I A stands for the indicator function of the set A.
(ii)-(a) Taking into account (i) and the fact that the p.g.f. of the variable ϕ 0 (k) is h k (u) = (1 − ρ + ρu) k , we have that there exists a positive constant C 1 such that Combining this with the expression of the conditional expectation we have that, there exist positive constants C 2 and C 3 such that In order to prove (ii)-(b), we can choose 0 < ζ < λµρ − 1, so that 1 < η = λµρ − ζ. Let us also define the sets A n = {ηY n < Y n+1 }, n ∈ N 0 . Then, it is not difficult to verify that, Moreover, for each l ∈ N 0 fixed, let us define the sets and now, we shall obtain a lower bound for P (A 0 |Y 0 = i). Note that since, f (1 − ρ + ρλ) i → 0, as i → ∞, then there exist I 0 ∈ N and 0 < ǫ < ζ such that ρλf As a consequence, on the one hand, by using Chebyshev's inequality, for i ≥ I 0 , On the other hand, using (a) we have that there exist positive constants K 1 and K 2 satisfying Note that since η > 1, by taking y ≥ I 0 we have yη l ≥ I 0 , for all l ∈ N 0 . Thus, from all the above we obtain Finally, we observe that N is a class of communicating states for this process and hence, using the same arguments as in [19], p.47, we can prove that P (Y n → ∞|Y 0 = y) > 0 for each y ∈ N.

B Proofs of the results in Section 4
In the following results we make use of the next lemma. The proof is easy and it is omitted.

Proof of Proposition 4.3
We only give the proofs for the predator females. The proofs are similar for the remaining conditional moments.
(i) Using the fact that G n ⊆ H n and Lemma B.1 (ii) Taking into account that G n ⊆ H n and Lemma B.1, we have that = α 2 σ 2 (T n ,T n ) + α(1 − α)ε(T n ,T n ) a.s.

Proof of Proposition 4.4
We only provide the proof for the number of females in the predator population; the remaining ones are similar.
(i) Taking into account F n ⊆ G n+1 , Proposition 4.3 (i), the expectation of the control variables, and the independence of T n andT n given Z n andZ n , we have that (ii) By Proposition 4.3 (i) and (ii) and bearing in mind that F n ⊆ G n+1 , we have that On the one hand, for the second and third terms by using Proposition 4.2 Analogously, for the first term Proof of Proposition 4.5 By condition (A1), Proposition 4.4 and the properties of the survival probability functions s(·) ands(·), we have that Analogously, we obtain that

Proof of Proposition 4.6
We provide the proofs for the number of predators; the arguments are similar for the number of preys.
In order to proof (ii) we observe that where A 1 = {Z n 0, Z n ∞}, and A 2 = {Z n 0,Z n ∞}. We shall start by proving that P (A 1 ) = 0. Given ω ∈ A 1 , then there exist 0 < B < ∞ such that for all n 0 , there exists n ≥ n 0 satisfying 0 < Z n (ω) ≤ B, i.e., then, and we conclude that P (A 1 ) = 0 by applying (i). Similar arguments lead to P (A 2 ) = 0, and thus, the result follows.

C Proofs of the results in Section 5
Proof of Proposition 5.1 First, observe that Now, taking into account the definition of the model and the fact that φ 0 (t, 0) = 0 a.s. for each t ∈ N 0 , by condition (A2) we conclude that for any n ∈ N 0 ,

Proof of Proposition 5.2
We follow the same arguments as in the proof of Lemma 2 in [17]. On the one hand, by Proposition 4.5 we have that E[Z n+1 |F n ] ≤αmZ n E[s(T n )|F n ] a.s. On the other hand, by the definition of the model and condition (A1), Z n = L(F n , M n ) ≤ F n ≤ T n , for each n ∈ N, and taking into account that the functions(·) is strictly decreasing we have that, for all n ≥ 1,s(Z n ) ≥s(T n ), and then Since lim x→∞s (x) = 0, there exists A > 0 such that |x| ≥ A impliess(x) ≤ 1 αm , and then E[Z n+1 |F n ] ≤Z n a.s. on {Z n ≥ A}.
We shall prove now that this implies that P (i,j) (Z n → ∞,Z n → ∞) = 0. To that end, we shall prove that, for every N > 0, Let us fix N > 0 and define the stopping time: if inf n≥N Z n ≥ A, min{n ≥ N : Z n < A}, otherwise, and define also the sequence of r.v.s {Y n } n∈N 0 as follows: To obtain (8), we show that {Y n } n∈N 0 is a non-negative super-martingale with respect to {F N +n } n∈N 0 . Indeed, the variable Y n is F N +n -measurable for any n ≥ 0.
Let us fix n ≥ 0, if Z N +k ≥ A, for each k = 0, . . . , n, then T (A) ≥ N + n + 1 and from (7) we obtain that Applying the martingale convergence theorem, we obtain the almost sure convergence of the sequence {Y n } n∈N 0 to a non-negative and finite limit, Y ∞ where Y ∞ = lim n→∞Zn , if inf n≥N Z n > A, Z T (A) , otherwise.
Therefore we deduce (8) and the proof is completed.

D Proofs of the results in Section 6
Proof of Proposition 6.2 From Proposition 4.5 it is immediate to see that if ρ 2 mα ≤ 1 andρ 2mα ≤ 1, then {Z n } n∈N 0 and {Z n } n∈N 0 are both non-negative supermartingales with respect to the σ-algebra F n and therefore both of them converge a.s. to a finite limit. Now, the result is derived by Proposition 4.7 (ii).
The result follows with the same arguments as in the proof of Proposition 5.2.
(ii) First, notice that {T n → ∞,T n → ∞} = {Z n → ∞,Z n → ∞} a.s. and thus, it is enough to prove that P (i,j) (T n → ∞,T n → ∞) > 0 if min{ρ 2 mα,ρ 1mα } > 1. To do that, first we prove that there exist I 1 , J 1 ∈ N 0 large enough so that P (i,j) (T n → ∞,T n → ∞) > 0, for each i, j ∈ N satisfying i ≥ I 1 and j ≥ J 1 , and then we prove that this implies the desired result.
Then, it is not difficult to verify that, P (i,j) (T n → ∞,T n → ∞) ≥ P (i,j) (∩ ∞ n=0 A n ) = lim n→∞ P (i,j) (A 0 ) n l=1 P A l | ∩ l−1 k=0 A k ∩ {Z 0 = i,Z 0 = j} , and for each l ∈ N fixed we have that As a consequence, we only need to obtain convenient lower bounds for P (i,j) (A 0 ) and P (A 1 |(T 1 ,T 1 ) = (i 0 , j 0 )). For the former we have For the latter, we shall obtain a convenient upper bound for P (A c 1 |(T 1 ,T 1 ) = (i 0 , j 0 )). First, from Proposition 4.6 we have and then P (i,j) (Z n → ∞,Z n → ∞) > 0 for i ≥ I 1 and j ≥ J 1 . Finally, we prove that the result also holds for any i, j ∈ N. Let us fix i, j ∈ N such that either i < I 1 or j < J 1 . By Proposition 4.1 (ii), there exists k 0 ∈ N such that P (Z k 0 = I 1 ,Z k 0 = J 1 |Z 0 = i,Z 0 = j) > 0. Then, we have P (Z n → ∞,Z n → ∞|Z 0 = i,Z 0 = j) ≥ P (Z n → ∞,Z n → ∞|Z k 0 = I 1 ,Z k 0 = J 1 ) · P (Z k 0 = I 1 ,Z k 0 = J 1 |Z 0 = i,Z 0 = j) > 0.