A Predator–Prey Two-Sex Branching Process

Abstract

In this paper, we present the first stochastic process to describe the interaction of predator and prey populations with sexual reproduction. Specifically, we introduce a two-type two-sex controlled branching model. 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. In view of their potential for the conservation of species, we provide necessary and sufficient conditions for the ultimate extinction of both species, the fixation of one of them 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.

1. Introduction

Predator–prey models have been widely studied in the literature since the introduction of the first model in [1]. This process aimed at modelling the trophic interactions between both species and as later in [2], non-linear differential equations were used to describe the interaction of a predator–prey 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, [3,4,5,6,7], which deal with theoretical questions, or [8,9], 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, [10] 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 [11] 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 [11] is performed in [12]. In [13], 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 [14] and it focusses on the study of necessary and sufficient conditions for the fixation of both populations. Later, [15] analyses a counterpart of the predator–prey model in [13] in the continuous-time setting by making use of linear birth-death processes. More recently, in [16] 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 [17] where a branching random walk is used to describe predator–prey populations.

Despite the broad literature on this topic, minor attempts have been made to include sexual reproduction in predator–prey models. For instance, [18] proposed a simulated model for a very specific situation in a genetic context where predators have sexual reproduction, preys can be sexual or asexual, the progenies of both species follow binomial distributions and the predator and prey population sizes are fixed. However, this last assumption is not common in the nature and then, the main weakness of the mentioned model is the lack of realism in terms of population dynamics. Apart from this paper, to the best of our knowledge none of the previous stochastic models covers the fact that both species have a sexual reproduction. This 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, [19] studied populations of sea lions and penguins as predators and preys, respectively, and in [20], a social network of giraffe populations is studied bearing in mind the presence of lions which preferentially prey on giraffe calves. Apart from the aforementioned predator–prey populations, there is a large number of predator–prey systems that have recently drawn the attention of biologists because of their -direct or indirect- impact on human life. Some examples are the Atlantic Salmon and its predator, the European catfish (where the number of the former are declining, see [21]), the American blue crab that preys on prawns and that has spread throughout the Mediterranean affecting commercial fishing (see [22]) or the Chinook salmon, whose population size has decreased as a consequence of conservation strategies for resident killer whales (see [23]) among others. The reader is referred to [24,25,26,27,28], or [29] for further examples where our model could be fitted.

Our aim in this paper is to introduce a predator–prey stochastic process by using the two-sex branching processes theory (see [30]) 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. A major advantage of the present model versus other deterministic predator–prey models is that our model enables to describe the random nature of the reproduction, mating and interaction among individuals, as is usual in the nature. Furthermore, for the first time a stochastic predator–prey model includes the sexual reproduction among individuals of both species. This represents a clear progress in this line of research given that the percentage of eukaryote organisms with sexual reproduction is estimated in 99.99% approximately (see [31]). Due to the variety of mating systems, 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, [32,33,34,35,36] or [37]).

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 natural 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. As a starting point for the study of these processes and motivated by their importance in conservation of species, we examine conditions for one of the species -predator or prey- to become extinct or to have a positive probability of survival. These conditions might help scientists to design strategies for the conservation of endanger species in real settings. 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 seven 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 Section 5 and Section 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 in Section 2, Section 4, Section 5, and Section 6 are gathered in Appendix A, Appendix B, Appendix C, and Appendix D, respectively.

2. A Controlled Two-Sex Branching Process

Before providing the definition of the predator–prey two-sex branching process, in this section we introduce 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 [30] 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. Moreover, this model differs from the one introduced in [38] in the fact that in our case the control is not made on the number of couples, but on the total number of individuals, which seems to be a more reasonable assumption. First, we shall give an intuitive insight of the model.

Throughout this section, the variable $Y_n$ represents the number of couples at generation n, while $X_n$ is the total number of individuals—females and males—in the population at that generation and the function ${\phi }_n(·)$ denotes the control on these individuals. In this model, we distinguish three consecutive phases at each generation: 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 ${\phi }_n\left(X_n\right)$ at generation n. Each one of these individuals can be female or male. It is a female with probability $\lambda$ and it is a male with probability $1-\lambda$. We also write ${\overline{F}}_n$ and ${\overline{M}}_n$ to refer to the number of females and males, respectively, that survive the control phase at generation n. In particular, ${\phi }_n\left(X_n\right)={\overline{F}}_n+{\overline{M}}_n$. The second phase is the mating phase where ${\overline{F}}_n$ females and ${\overline{M}}_n$ males at generation n mate to produce $Y_n$ couples by some mating system L, $Y_n=L({\overline{F}}_n,{\overline{M}}_n)$. The last phase is a reproduction phase where the $Y_n$ couples give birth to their offspring. We denote the number of offspring given by the i-th couple in the n-th generation by ${\xi }_{ni}$, and consequently, all these offspring constitute the total number of individuals at the following generation, $X_{n+1}$. Next, we provide the formal definition of the model.

Let us assume that the families of random variables (r.v.s) $\{{\xi }_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$ and $\{{\phi }_n\left(k\right):n,k\in {\mathbb{N}}_0\}$ are independent. 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 $\{{\phi }_n\left(k\right):k\in {\mathbb{N}}_0\}$, $n\in {\mathbb{N}}_0$, are independent stochastic processes with the same one-dimensional distribution. Let us also consider the sequences ${\{Y_n\}}_{n\in {\mathbb{N}}_0}$ and ${\{X_n\}}_{n\in \mathbb{N}}$ defined as:

$$Y_0=y_0\in \mathbb{N},X_{n+1}=\sum _{i=1}^{Y_n}{\xi }_{ni},Y_{n+1}=L({\overline{F}}_{n+1},{\overline{M}}_{n+1})n\in {\mathbb{N}}_0,$$

(1)

where ${\mathbb{N}}_0=\mathbb{N}\cup \{0\}$, $L:{\mathbb{N}}_0^2\to {\mathbb{N}}_0$ is a deterministic mating function and $({\overline{F}}_{n+1},{\overline{M}}_{n+1})$ is a random vector that follows a multinomial distribution with parameters y and $(\lambda ,1-\lambda )$ conditionally on $\{{\phi }_{n+1}\left(X_{n+1}\right)=y\}$, with $0<\lambda <1$. The process ${\{Y_n\}}_{n\in {\mathbb{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:

(C1)

$L(x,y)\le x$, for each $x,y\in {\mathbb{N}}_0$.

(C2)

$L(0,y)=L(x,0)=0$, for each $x,y\in {\mathbb{N}}_0$.

Henceforth, we consider the promiscuous mating function as the mating function L, which is defined as $L(x,y)=xmin\{1,y\}$. We denote the probability distribution of the variable ${\xi }_{01}$ by ${\{q_k\}}_{k\in {\mathbb{N}}_0}$, which is named the offspring distribution or reproduction law, and to avoid trivialities it is assumed to satisfy $q_0+q_1+q_2<1$.

Remark 1.

(a) 

Condition (C1) is a natural hypothesis and it establishes that the number of couples is always less than or equal to the number of females. The meaning of the Condition (C2) is obvious, if there is no female or male in the population, no couple can be formed.

(b) 

Apart from the promiscuous mating other classic examples of mating functions used in two-sex branching processes theory and satisfying Conditions (C1) and (C2) are the following:

(i) 

Polygamous mating: in this situation females and males could form couples with more than one individual of the opposite sex. It can be modelled through the function $L(x,y)=min\{x,dy\},d>1$.

(ii) 

Perfect fidelity mating: in this case each individual mates with only one individual of the opposite sex forming an exclusive couple. The mating function is $L(x,y)=min\{x,y\}$.

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\in {\mathbb{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 ${\phi }_0\left(0\right)=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 [39] 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.

Proposition 1.

Let ${\{Y_n\}}_{n\in {\mathbb{N}}_0}$ be a BBPCI. Then:

(i) 

$\mathbb{N}$ is a family of transient states.

(ii) 

The classic duality extinction-explosion holds, that is, $P(Y_n\to 0)+P(Y_n\to \infty )=1$.

In order to establish the following results let us denote the mean and variance of the variable ${\xi }_{01}$ as $0<\mu <\infty$, and $0<{\delta }^2<\infty$, respectively, and write $ϵ\left(k\right)=E\left[{\phi }_0\left(k\right)\right]$ and ${\delta }^2\left(k\right)=\mathrm{Var}\left[{\phi }_0\left(k\right)\right]$, and assume that both of them are positive and finite for each $k\in {\mathbb{N}}_0$.

Theorem 1.

Let ${\{Y_n\}}_{n\in {\mathbb{N}}_0}$ be a BBPCI. If $\lambda \mu ϵ\left(k\right)\le k$, for each $k\in {\mathbb{N}}_0$, then $P(Y_n\to \infty |Y_0=y)=0$, for each $y\in \mathbb{N}$.

An immediate consequence of the previous theorem is that a controlled two-sex branching process with control functions ${\phi }_0\left(k\right)$ following binomial distributions with parameters k and $0<\rho <1$ becomes extinct a.s. if $\lambda \mu \rho \le 1$. This remark will be useful in the proof of the results in Section 6.

Theorem 2.

Let ${\{Y_n\}}_{n\in {\mathbb{N}}_0}$ be a BBPCI.

(i) 

For each $n\in {\mathbb{N}}_0$, the conditional expectation is:

$$E[Y_{n+1}|Y_n]=\lambda E[ϵ\left(X_{n+1}\right)|Y_n]-\lambda E[h_{X_{n+1}}^{\prime }\left(\lambda \right)|Y_n]a.s.,$$

where $h_k(·)$ denotes the probability generating function (p.g.f.) of ${\phi }_0\left(k\right)$.

(ii) 

If the control function ${\phi }_0\left(k\right)$ follows a binomial distribution with parameters k and ρ, then:

(a) 

If $f(·)$ denotes the p.g.f of the variable ${\xi }_{01}$, then there exist constants $C_1,C_2,C_3>0$ such that for $n\in {\mathbb{N}}_0$,

$$\begin{array}{cc}\hfill E[Y_{n+1}|Y_n]& =\lambda \rho \mu Y_n-C_1{Y}_{n}f{(1-\rho +\rho \lambda )}^{Y_n}a.s.,\hfill \\ \hfill \mathrm{Var}[Y_{n+1}|Y_n]& \le C_2{Y}_n+C_3{Y}_n^{2}f{(1-\rho +\rho \lambda )}^{Y_n}a.s.\hfill \end{array}$$

(b) 

If $\lambda \mu \rho >1$, then $P(Y_n\to \infty |Y_0=y)>0$, for each $y\in \mathbb{N}$.

3. Definition of a Predator–Prey Two-Sex Branching Model

Having described the controlled two-sex branching process, in this section we introduce a predator–prey 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 and then we provide its intuitive description.

Let $\{t_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$, $\{{\tilde{t}}_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$, $\{{\varphi }_n(t,\tilde{t}):n,t,\tilde{t}\in {\mathbb{N}}_0\}$ and $\{{\tilde{\varphi }}_n(t,\tilde{t}):n,t,\tilde{t}\in {\mathbb{N}}_0\}$ be independent families of non-negative and integer valued r.v.s defined on the same probability space $(\mathrm{\Omega },\mathcal{A},P)$, and assume that:

(i)

The r.v.s of the sequence $\{t_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$ are i.i.d. with probability distribution $p={\{p_k\}}_{k\in {\mathbb{N}}_0}$, where $p_k=P(t_{01}=k)$, $k\in {\mathbb{N}}_0$.

(ii)

The r.v.s of the family $\{{\tilde{t}}_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$ are i.i.d. with probability distribution $\tilde{p}={\{{\tilde{p}}_k\}}_{k\in {\mathbb{N}}_0}$, where ${\tilde{p}}_k=P({\tilde{t}}_{01}=k)$, $k\in {\mathbb{N}}_0$.

(iii)

For each $n,t,\tilde{t}\in {\mathbb{N}}_0$, ${\varphi }_n(t,\tilde{t})$ is a r.v. following a binomial distribution with parameters t and $s\left(\tilde{t}\right)$, where $s:\mathbb{R}\to [0,1]$ is a strictly increasing function which is continuous at 0 and such that, for certain $0\le {\rho }_1<{\rho }_2\le 1$,

$$s\left(0\right)={\rho }_1,\mathrm{and}\underset{\tilde{t}\to \infty }{lim}s\left(\tilde{t}\right)={\rho }_2.$$

(2)

(iv)

For each $n,t,\tilde{t}\in {\mathbb{N}}_0$, ${\tilde{\varphi }}_n(t,\tilde{t})$ is a r.v. following a binomial distribution with parameters $\tilde{t}$ and $\tilde{s}\left(t\right)$, where $\tilde{s}:\mathbb{R}\to [0,1]$ is a strictly decreasing function which is continuous at 0 and such that, for certain $0\le {\tilde{\rho }}_1<{\tilde{\rho }}_2\le 1$,

$$\tilde{s}\left(0\right)={\tilde{\rho }}_2,\mathrm{and}\underset{t\to \infty }{lim}\tilde{s}\left(t\right)={\tilde{\rho }}_1.$$

(3)

A predator–prey two-sex branching process is a bivariate stochastic process ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ defined recursively as:

$$(Z_0,{\tilde{Z}}_0)=(N,\tilde{N})\in {\mathbb{N}}^2,(T_{n+1},{\tilde{T}}_{n+1})=\left(\sum _{i=1}^{Z_n}{t}_{ni},\sum _{i=1}^{{\tilde{Z}}_n}{\tilde{t}}_{ni}\right),n\in {\mathbb{N}}_0,$$

(4)

and

$$(Z_{n+1},{\tilde{Z}}_{n+1})=\left(L(F_{n+1},M_{n+1}),\tilde{L}({\tilde{F}}_{n+1},{\tilde{M}}_{n+1})\right),n\in {\mathbb{N}}_0,$$

where conditionally on $\{{\varphi }_{n+1}(T_{n+1},{\tilde{T}}_{n+1})=k\}$, the random vector $(F_{n+1},M_{n+1})$ follows a multinomial distribution with parameters k and $(\alpha ,1-\alpha )$, with $0<\alpha <1$, and conditionally on $\{{\tilde{\varphi }}_{n+1}(T_{n+1},{\tilde{T}}_{n+1})=\tilde{k}\}$ the random vector $({\tilde{F}}_{n+1},{\tilde{M}}_{n+1})$ follows a multinomial distribution with parameters $\tilde{k}$ and $(\tilde{\alpha },1-\tilde{\alpha })$, with $0<\tilde{\alpha }<1$. Moreover, the empty sums in (4) are assumed to be 0.

The process defined above enables us to model a predator–prey system with non-overlap** 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$ and ${\tilde{T}}_n$ represent the total number of predators and preys, respectively, at generation n, whereas $F_n$ and $M_n$ (${\tilde{F}}_n$ and ${\tilde{M}}_n$) are the total numbers of progenitor predator females and predator males (prey females and prey males), respectively, at generation n. Moreover, $Z_n$ and ${\tilde{Z}}_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 r.v.s $\{t_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$, and $\{{\tilde{t}}_{ni}:n\in {\mathbb{N}}_0,i\in \mathbb{N}\}$, where $t_{ni}$ denotes the number of offspring of the i-th predator couple at generation n while ${\tilde{t}}_{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\in {\mathbb{N}}_0}$ and $\tilde{p}={\{{\tilde{p}}_k\}}_{k\in {\mathbb{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 ${\sigma }^2$, respectively, for the predators and $\tilde{m}$ and ${\tilde{\sigma }}^2$ for the preys. We recall that we consider that the mating functions L and $\tilde{L}$ are the promiscuous mating function and hence, in order to avoid trivialities, we also assume that $p_0+p_1+p_2<1$ and ${\tilde{p}}_0+{\tilde{p}}_1+{\tilde{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, ${\tilde{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 and ${\tilde{T}}_{n+1}=\tilde{t}$ preys in the population, the number of predators and preys which survive are given by the r.v.s ${\varphi }_{n+1}(t,\tilde{t})$ and ${\tilde{\varphi }}_{n+1}(t,\tilde{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 ${\varphi }_{n+1}(t,\tilde{t})$ and ${\tilde{\varphi }}_{n+1}(t,\tilde{t})$ are binomial distributions.

More specifically, ${\varphi }_{n+1}(t,\tilde{t})$ is assumed to follow a binomial distribution with parameters t and $s\left(\tilde{t}\right)$, where the $s\left(\tilde{t}\right)$ represents the survival probability of a predator given that there are $\tilde{t}$ preys in the population. The condition on the monotonicity of the function $s:\mathbb{R}\to [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 ${\rho }_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, ${\rho }_1=0$ implies the extinction of the predator population when there are not any preys in the population at some generation. Moreover, ${\rho }_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 ${\rho }_2=1$ indicates the survival of all the predators when the number of preys goes to infinity.

Analogously, we assume that ${\tilde{\varphi }}_{n+1}(t,\tilde{t})$ follows a binomial distribution with parameters $\tilde{t}$ and $\tilde{s}\left(t\right)$, with $\tilde{s}\left(t\right)$ representing the survival probability of a prey given that there are t predators in the population. The condition on the monotonicity of the function $\tilde{s}:\mathbb{R}\to [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 ${\tilde{\rho }}_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 ${\tilde{\rho }}_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, ${\tilde{\rho }}_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 ${\tilde{\rho }}_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 $\alpha$ 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 $(\alpha ,1-\alpha )$, given that $T_{n+1}=t$, ${\tilde{T}}_{n+1}=\tilde{t}$ and ${\varphi }_{n+1}(t,\tilde{t})=y$. Similarly, there are ${\tilde{F}}_{n+1}$ females and ${\tilde{M}}_{n+1}$ males within the survivor prey population at generation $n+1$, and consequently, $({\tilde{F}}_{n+1},{\tilde{M}}_{n+1})$ follows a multinomial distribution of parameters $\tilde{y}$, and $(\tilde{\alpha },1-\tilde{\alpha })$, given that $T_{n+1}=t$, ${\tilde{T}}_{n+1}=\tilde{t}$ and ${\tilde{\varphi }}_{n+1}(t,\tilde{t})=\tilde{y}$, and where $\tilde{\alpha }$ 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}$ and ${\tilde{Z}}_{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 2.

(i) 

One can propose several functions $s(·)$ and $\tilde{s}(·)$ satisfying the previous assumptions. For example, for the survival of the predator population

$$s\left(x\right)={\rho }_2\left(1-a^{-x}\right)+{\rho }_1{a}^{-x},a>1,s\left(x\right)={\rho }_2\frac{x^k+{\rho }_1/{\rho }_2}{x^k+1},k>0,$$

and similarly, for the survival of the prey population

$$\tilde{s}\left(x\right)=\left({\tilde{\rho }}_2-{\tilde{\rho }}_1\right)a^{-x}+{\tilde{\rho }}_1,a>1,\tilde{s}\left(x\right)=1-\left(1-{\tilde{\rho }}_1\right)\frac{x^k+\left(1-{\tilde{\rho }}_2\right)/\left(1-{\tilde{\rho }}_1\right)}{x^k+1},k>0.$$

(ii) 

The presence of the mating phase makes difficult to compare the present stochastic model with its deterministic counterpart without including the sexual reproduction -the Lotka–Volterra (predator–prey) equations-. However, we can analyse the relation between them by examining the deterministic evolution of the number of individuals of each species that survive in presence of the other one. We note that $xs\left(\tilde{x}\right)$ is the expected number of predator survivors after the interaction with preys given that there are x predators and $\tilde{x}$ preys in the ecosystem. This mean number of predators corresponds to the interaction term in the Lotka–Volterra equations, while the term corresponding to the death of predators by natural causes is modelled via the reproduction phase. A similar comparison can be done for the preys. The mean number of prey survivors after the interaction with predators given that there are x predators and $\tilde{x}$ preys in the ecosystem is $\tilde{x}\tilde{s}\left(x\right)=\tilde{x}-\tilde{x}\tilde{d}\left(x\right)$, where $\tilde{s}(·)=1-\tilde{d}(·)$, and $\tilde{d}(·)$ represents the probability that a prey dies as a consequence of the presence of predators. Thus, the first term above corresponds to the growth of the number of preys in the Lotka–Volterra equations, while the second term—expected number of preys that die because of the predators—corresponds to the interaction term in the Lotka–Volterra equations. We also remark that the increment in the number of preys due to the reproduction is modelled through the reproduction phase as described above.

4. 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,{\tilde{Z}}_n)\}}_{n\in {\mathbb{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 (C2).

In the following easy-to-prove proposition we state some properties of the states associated with the process.

Proposition 2.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then,

(i) 

Every non-null state $(i,j)\ne (0,0)$ is transient.

(ii) 

If $p_0+p_1+p_2<1$, ${\tilde{p}}_0+{\tilde{p}}_1+{\tilde{p}}_2<1$, and $0<{\rho }_1<{\rho }_2<1$, then the sets $\{(i,0):i>0\}$, $\{(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 process can move from the last set to the others in one step.

(iii) 

If $p_0+p_1+p_2<1$, ${\tilde{p}}_0+{\tilde{p}}_1+{\tilde{p}}_2<1$ and $0={\rho }_1<{\rho }_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,\tilde{x}\in {\mathbb{N}}_0$, the conditional expectations of the control variables are

$$\begin{array}{cc}\hfill \epsilon (x,\tilde{x})& :=E[{\varphi }_n(T_n,{\tilde{T}}_n)|T_n=x,{\tilde{T}}_n=\tilde{x}]=xs\left(\tilde{x}\right),\hfill \\ \hfill \tilde{\epsilon }(x,\tilde{x})& :=E[{\tilde{\varphi }}_n(T_n,{\tilde{T}}_n)|T_n=x,{\tilde{T}}_n=\tilde{x}]=\tilde{x}\tilde{s}\left(x\right).\hfill \end{array}$$

and the conditional variances of the control variables are

$$\begin{array}{cc}\hfill {\sigma }^2(x,\tilde{x})& :=Var[{\varphi }_n(T_n,{\tilde{T}}_n)|T_n=x,{\tilde{T}}_n=\tilde{x}]=xs\left(\tilde{x}\right)(1-s\left(\tilde{x}\right)),\hfill \\ \hfill {\tilde{\sigma }}^2(x,\tilde{x})& :=Var[{\tilde{\varphi }}_n(T_n,{\tilde{T}}_n)|T_n=x,{\tilde{T}}_n=\tilde{x}]=\tilde{x}\tilde{s}\left(x\right)(1-\tilde{s}\left(x\right)).\hfill \end{array}$$

For the next results, we introduce the following notation concerning the $\sigma$-algebras generated by the variables involved in the definition of the model. In particular, we denote

$$\begin{array}{cc}\hfill {\mathcal{F}}_n& =\sigma (Z_l,{\tilde{Z}}_l:l=0,\dots ,n),n\in {\mathbb{N}}_0,\hfill \\ \hfill {\mathcal{G}}_n& =\sigma (Z_l,{\tilde{Z}}_l,T_{l+1},{\tilde{T}}_{l+1}:l=0,\dots ,n-1),n\in \mathbb{N},\hfill \\ \hfill {\mathcal{H}}_n& =\sigma (Z_l,{\tilde{Z}}_l,T_{l+1},{\tilde{T}}_{l+1},{\varphi }_{l+1}(T_{l+1},{\tilde{T}}_{l+1}),{\tilde{\varphi }}_{l+1}(T_{l+1},{\tilde{T}}_{l+1}):l=0,\dots ,n-1),n\in \mathbb{N}.\hfill \end{array}$$

Then, we have that ${\mathcal{F}}_{n-1}\subseteq {\mathcal{G}}_n\subseteq {\mathcal{H}}_n$, for $n\in \mathbb{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.

Proposition 3 (Conditional moments of the number of individuals).

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for each $n\in {\mathbb{N}}_0$,

(i) 

$E[T_{n+1}|{\mathcal{F}}_n]=m{Z}_n$, and $E[{\tilde{T}}_{n+1}|{\mathcal{F}}_n]=m{\tilde{Z}}_n$.

(ii) 

$Var[T_{n+1}|{\mathcal{F}}_n]={\sigma }^2{Z}_n$, and $Var[{\tilde{T}}_{n+1}|{\mathcal{F}}_n]={\tilde{\sigma }}^2{\tilde{Z}}_n$.

Next, we establish some results concerning to the conditional moments of the total number of female and male predators and female and male preys.

Proposition 4 (Conditional moments of the number of females and males).

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for each $n\in {\mathbb{N}}_0$,

(i) 

The expected number of predator and prey females given the number of individuals and couples are

$$E[F_n|{\mathcal{G}}_n]=\alpha T_{n}s\left({\tilde{T}}_n\right),andE[{\tilde{F}}_n|{\mathcal{G}}_n]=\tilde{\alpha }{\tilde{T}}_n\tilde{s}\left(T_n\right),$$

and the corresponding to the number of predator and prey males are

$$E[M_n|{\mathcal{G}}_n]=(1-\alpha )T_{n}s\left({\tilde{T}}_n\right),andE[{\tilde{M}}_n|{\mathcal{G}}_n]=(1-\tilde{\alpha }){\tilde{T}}_n\tilde{s}\left(T_n\right).$$

(ii) 

The conditional variances of the number of predator and prey females and males given the number of individuals and couples are

$$\begin{array}{cc}\hfill Var[F_n|{\mathcal{G}}_n]& ={\alpha }^2{T}_{n}s\left({\tilde{T}}_n\right)(1-s\left({\tilde{T}}_n\right))+\alpha (1-\alpha )T_{n}s\left({\tilde{T}}_n\right),\hfill \\ \hfill Var[{\tilde{F}}_n|{\mathcal{G}}_n]& ={\tilde{\alpha }}^2{\tilde{T}}_n\tilde{s}\left(T_n\right)(1-\tilde{s}\left(T_n\right))+\tilde{\alpha }(1-\tilde{\alpha }){\tilde{T}}_n\tilde{s}\left(T_n\right),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill Var[M_n|{\mathcal{G}}_n]& ={(1-\alpha )}^2{T}_{n}s\left({\tilde{T}}_n\right)(1-s\left({\tilde{T}}_n\right))+\alpha (1-\alpha )T_{n}s\left({\tilde{T}}_n\right),\hfill \\ \hfill Var[{\tilde{M}}_n|{\mathcal{G}}_n]& ={(1-\tilde{\alpha })}^2{\tilde{T}}_n\tilde{s}\left(T_n\right)(1-\tilde{s}\left(T_n\right))+\tilde{\alpha }(1-\tilde{\alpha }){\tilde{T}}_n\tilde{s}\left(T_n\right),\hfill \end{array}$$

respectively.

Proposition 5 (Conditional moments of the number of females and males).

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for each $n\in {\mathbb{N}}_0$,

(i) 

The expected number of predator and prey females given the number of couples are

$$E[F_{n+1}|{\mathcal{F}}_n]=\alpha m{Z}_{n}E[s\left({\tilde{T}}_{n+1}\right)|{\mathcal{F}}_n],andE[{\tilde{F}}_{n+1}|{\mathcal{F}}_n]=\tilde{\alpha }\tilde{m}{\tilde{Z}}_{n}E[\tilde{s}\left(T_{n+1}\right)|{\mathcal{F}}_n],$$

and the corresponding to the number of predator and prey males are

$$E[M_{n+1}|{\mathcal{F}}_n]=(1-\alpha )m{Z}_{n}E[s\left({\tilde{T}}_{n+1}\right)|{\mathcal{F}}_n],andE[{\tilde{M}}_{n+1}|{\mathcal{F}}_n]=(1-\tilde{\alpha })\tilde{m}{\tilde{Z}}_{n}E[\tilde{s}\left(T_{n+1}\right)|{\mathcal{F}}_n].$$

(ii) 

The conditional variance of the number of predator and prey females and males given the number of couples are

$$\begin{array}{cc}\hfill Var[F_{n+1}|{\mathcal{F}}_n]& ={\alpha }^2{m}^2{Z}_n^{2}Var[s\left({\tilde{T}}_{n+1}\right)|{\mathcal{F}}_n]+\alpha m{Z}_{n}E[s\left({\tilde{T}}_{n+1}\right)|{\mathcal{F}}_n]+{\alpha }^2{Z}_n({\sigma }^2-m)E[s{\left({\tilde{T}}_{n+1}\right)}^2|{\mathcal{F}}_n],\hfill \\ \hfill Var[{\tilde{F}}_{n+1}|{\mathcal{F}}_n]& ={\tilde{\alpha }}^2{\tilde{m}}^2{\tilde{Z}}_n^{2}Var[\tilde{s}\left(T_{n+1}\right)|{\mathcal{F}}_n]+\tilde{\alpha }\tilde{m}{\tilde{Z}}_{n}E[\tilde{s}\left(T_{n+1}\right)|{\mathcal{F}}_n]+{\tilde{\alpha }}^2{\tilde{Z}}_n({\tilde{\sigma }}^2-\tilde{m})E[\tilde{s}{\left(T_{n+1}\right)}^2|{\mathcal{F}}_n].\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill Var[M_{n+1}|{\mathcal{F}}_n]& ={(1-\alpha )}^2{m}^2{Z}_n^{2}Var[s\left({\tilde{T}}_{n+1}\right)|{\mathcal{F}}_n]+(1-\alpha )m{Z}_{n}E[s\left({\tilde{T}}_{n+1}\right)|{\mathcal{F}}_n]\hfill \\ \hfill & \phantom{=}+{(1-\alpha )}^2{Z}_n({\sigma }^2-m)E[s{\left({\tilde{T}}_{n+1}\right)}^2|{\mathcal{F}}_n],\hfill \\ \hfill Var[{\tilde{M}}_{n+1}|{\mathcal{F}}_n]& ={(1-\tilde{\alpha })}^2{\tilde{m}}^2{\tilde{Z}}_n^{2}Var[\tilde{s}\left(T_{n+1}\right)|{\mathcal{F}}_n]+(1-\tilde{\alpha })\tilde{m}{\tilde{Z}}_{n}E[\tilde{s}\left(T_{n+1}\right)|{\mathcal{F}}_n]\hfill \\ \hfill & \phantom{=}+{(1-\tilde{\alpha })}^2{\tilde{Z}}_n({\tilde{\sigma }}^2-\tilde{m})E[\tilde{s}{\left(T_{n+1}\right)}^2|{\mathcal{F}}_n],\hfill \end{array}$$

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.

Proposition 6.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for each $n\in {\mathbb{N}}_0$,

$$E[Z_{n+1}|{\mathcal{F}}_n]\le \alpha m{Z}_{n}E[s\left({\tilde{T}}_n\right)|{\mathcal{F}}_n]\le {\rho }_2\alpha m{Z}_n,andE[{\tilde{Z}}_{n+1}|{\mathcal{F}}_n]\le \tilde{\alpha }\tilde{m}{\tilde{Z}}_{n}E[\tilde{s}\left(T_n\right)|{\mathcal{F}}_n]\le {\tilde{\rho }}_2\tilde{\alpha }\tilde{m}{\tilde{Z}}_n.$$

Remark 3.

Since ${\rho }_1\le s\left(t\right)\le {\rho }_2$, and ${\tilde{\rho }}_1\le \tilde{s}\left(t\right)\le {\tilde{\rho }}_2$, for all $t\ge 0$, we obtain upper and lower bounds for the previous conditional expectations,

$$\begin{array}{ccc}\hfill {\rho }_1\alpha T_n\le & E[F_n|{\mathcal{G}}_n]& \le {\rho }_2\alpha T_n,\hfill \\ \hfill {\rho }_1\alpha m{Z}_n\le & E[F_{n+1}|{\mathcal{F}}_n]& \le {\rho }_2\alpha m{Z}_n.\hfill \end{array}$$

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.

Proposition 7 (Conditional moments of the number of individuals).

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for any $n\in \mathbb{N}$,

(i) 

$\begin{array}{c}E[T_{n+1}|{\mathcal{G}}_n]=\alpha m{T}_{n}s\left({\tilde{T}}_n\right)-\alpha m{T}_{n}s\left({\tilde{T}}_n\right){(1-s\left({\tilde{T}}_n\right)+\alpha s\left({\tilde{T}}_n\right))}^{T_n-1}.\hfill \\ E[{\tilde{T}}_{n+1}|{\mathcal{G}}_n]=\tilde{\alpha }\tilde{m}{\tilde{T}}_n\tilde{s}\left(T_n\right)-\tilde{\alpha }\tilde{m}{\tilde{T}}_n\tilde{s}\left(T_n\right){(1-\tilde{s}\left(T_n\right)+\tilde{\alpha }\tilde{s}\left(T_n\right))}^{{\tilde{T}}_n-1}.\hfill \end{array}$

(ii) 

$\begin{array}{c}Var[T_{n+1}|{\mathcal{G}}_n]\le ({\sigma }^2+m^2)\alpha T_{n}s\left({\tilde{T}}_n\right)+2{\alpha }^2{m}^2{T}_n^{2}s{\left({\tilde{T}}_n\right)}^2{(1-s\left({\tilde{T}}_n\right)+\alpha s\left({\tilde{T}}_n\right))}^{T_n-1}.\hfill \\ Var[{\tilde{T}}_{n+1}|{\mathcal{G}}_n]\le ({\tilde{\sigma }}^2+{\tilde{m}}^2)\tilde{\alpha }{\tilde{T}}_n\tilde{s}\left(T_n\right)+2{\tilde{\alpha }}^2{\tilde{m}}^2{\tilde{T}}_n^2\tilde{s}{\left(T_n\right)}^2{(1-\tilde{s}\left(T_n\right)+\tilde{\alpha }\tilde{s}\left(T_n\right))}^{{\tilde{T}}_n-1}.\hfill \end{array}$

To conclude this section we establish the following result where we derive the usual property of branching processes known as the extinction-explosion dichotomy.

Proposition 8.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then,

(i) 

$P({lim\; inf}_{n\to \infty }(Z_n,{\tilde{Z}}_n)=(k,\tilde{k}))=0$, and $P({lim\; sup}_{n\to \infty }(Z_n,{\tilde{Z}}_n)=(k,\tilde{k}))=0$, for each $(k,\tilde{k})\in {\mathbb{N}}_0^2\backslash \{(0,0)\}$.

(ii) 

$P(Z_n\to 0,{\tilde{Z}}_n\to 0)+P(Z_n\to \infty ,{\tilde{Z}}_n\to \infty )+P(Z_n\to \infty ,{\tilde{Z}}_n\to 0)+P(Z_n\to 0,{\tilde{Z}}_n\to \infty )=1$.

The sets $\{Z_n\to 0,{\tilde{Z}}_n\to 0\}$, $\{Z_n\to \infty ,{\tilde{Z}}_n\to \infty \}$, $\{Z_n\to \infty ,{\tilde{Z}}_n\to 0\}$ and $\{Z_n\to 0,{\tilde{Z}}_n\to \infty \}$ 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\to 0\}$, and $\{Z_n\to \infty \}$, and the extinction and survival of the prey population as $\{{\tilde{Z}}_n\to 0\}$, and $\{{\tilde{Z}}_n\to \infty \}$, then in view of Proposition 8 it is immediate that,

$$\begin{array}{cc}\hfill \{Z_n\to 0\}& =\{Z_n\to 0,{\tilde{Z}}_n\to 0\}\cup \{Z_n\to 0,{\tilde{Z}}_n\to \infty \}a.s.\hfill \\ \hfill \{Z_n\to \infty \}& =\{Z_n\to \infty ,{\tilde{Z}}_n\to 0\}\cup \{Z_n\to \infty ,{\tilde{Z}}_n\to \infty \}a.s.\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \{{\tilde{Z}}_n\to 0\}& =\{Z_n\to 0,{\tilde{Z}}_n\to 0\}\cup \{Z_n\to \infty ,{\tilde{Z}}_n\to 0\}a.s.\hfill \\ \hfill \{{\tilde{Z}}_n\to \infty \}& =\{Z_n\to 0,{\tilde{Z}}_n\to \infty \}\cup \{Z_n\to \infty ,{\tilde{Z}}_n\to \infty \}a.s.\hfill \end{array}$$

(6)

Remark 4.

Note that, in the cases of fixation of the predator and the prey, the process behaves as a BBPCI defined in (1) (see Section 2).

5. 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 [40]), Eastern Island (see [41]), Macquarie Island (see [42]), isolated regions in Finland and Northwest Russia (see [43]) or seafloor plateau (see [44]).

In terms of our model, this situation can be expressed as ${\rho }_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 ${\tilde{\rho }}_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 ${\rho }_2\le 1$ and ${\tilde{\rho }}_2\le 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, if ${\tilde{Z}}_n=0$ for some $n>0$, then ${\tilde{T}}_k=0$ and ${\tilde{Z}}_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\left(0\right)=0$, and consequently, ${\varphi }_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 since $\tilde{s}\left(0\right)={\rho }_2$, ${\tilde{\varphi }}_k(\tilde{t},0)$ follows a binomial distribution with parameters $\tilde{t}$ and ${\rho }_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,{\tilde{Z}}_0)=(i,j))$, in the remaining of this paper.

5.1. 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.

Proposition 9.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for any initial values $i,j>0$, $P_{(i,j)}(Z_n\to \infty ,{\tilde{Z}}_n\to 0)=0$.

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.

Proposition 10.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. Then, for any initial values $i,j>0$, $P_{(i,j)}(Z_n\to \infty ,{\tilde{Z}}_n\to \infty )=0$.

Note that by (5), the previous results and Proposition 8 (ii) we obtain the certain extinction of the predator population, that is, $P_{(i,j)}(Z_n\to 0)=1$, for any initial values $i,j>0$.

5.2. The Fixation of the Prey Population

From Propositions 9 and 10, 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 ${\tilde{\rho }}_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 [30] can be applied in this setting and we immediately obtain the following result.

Proposition 11.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. If ${\tilde{\rho }}_2=1$, for any initial values $i,j>0$, $P(Z_n\to 0,{\tilde{Z}}_n\to \infty )>0$ if and only if $\tilde{\alpha }\tilde{m}>1$.

On the other hand, if ${\tilde{\rho }}_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 ${\tilde{\rho }}_2$. In this case, the result is a direct consequence of Theorems 1 and 2 (ii) (b) provided in the mentioned section.

Proposition 12.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. If ${\tilde{\rho }}_2<1$, for any initial values $i,j>0$, $P_{(i,j)}(Z_n\to 0,{\tilde{Z}}_n\to \infty )>0$ if and only if $\tilde{\alpha }\tilde{m}{\tilde{\rho }}_2>1$.

6. Predator–Prey Non-Isolated System

In this section we consider the case that $0<{\rho }_1<{\rho }_2\le 1$ and $0<{\tilde{\rho }}_1<{\tilde{\rho }}_2<1$. Thus, contrary to Section 5, the predators have a positive probability of survival in absence of the prey population (${\rho }_1>0$) due to the availability of other food resources and also a limited appetite (${\tilde{\rho }}_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 (${\tilde{\rho }}_2<1$).

6.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 $\gamma$, where $\gamma ={\rho }_1$ in the case of the predator fixation and $\gamma ={\tilde{\rho }}_2$ in the case of the prey fixation. Thus, by using Theorems 1 and 2 (ii) (b) we have the following result.

Proposition 13.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. For any initial values $i,j>0$:

(i) 

$P_{(i,j)}(Z_n\to \infty ,{\tilde{Z}}_n\to 0)>0$ if and only if ${\rho }_1\alpha m>1$.

(ii) 

$P_{(i,j)}(Z_n\to 0,{\tilde{Z}}_n\to \infty )>0$ if and only if ${\tilde{\rho }}_2\tilde{\alpha }\tilde{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 13 is the following corollary:

Corollary 1.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. For any initial values $i,j>0$:

(i) 

If ${\rho }_1\alpha m>1$, then $P_{(i,j)}(Z_n\to \infty )>0$.

(ii) 

If ${\tilde{\rho }}_2\tilde{\alpha }\tilde{m}>1$, then $P_{(i,j)}({\tilde{Z}}_n\to \infty )>0$.

6.2. 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 14.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. For any initial values $i,j>0$:

(i) 

If ${\rho }_2\alpha m\le 1$, then $P_{(i,j)}(Z_n\to 0)=1$.

(ii) 

If ${\tilde{\rho }}_2\tilde{\alpha }\tilde{m}\le 1$, then $P_{(i,j)}({\tilde{Z}}_n\to 0)=1$.

Taking into account (5) and (6) and from Proposition 14, we deduce the following result on the coexistence of the species:

Corollary 2.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. For any initial values $i,j>0$, if $min\{{\rho }_2\alpha m,{\tilde{\rho }}_2\tilde{\alpha }\tilde{m}\}\le 1$, then $P_{(i,j)}(Z_n\to \infty ,{\tilde{Z}}_n\to \infty )=0$.

We note that there is always a positive probability for the complete extinction of the predator–prey 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$ and ${\tilde{p}}_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 15.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. For any initial values $i,j>0$, $P_{(i,j)}(Z_n\to 0,{\tilde{Z}}_n\to 0)=1$ if and only if $max\{{\rho }_{1}m\alpha ,{\tilde{\rho }}_2\tilde{m}\tilde{\alpha }\}\le 1$.

6.3. The Predator and Prey Coexistence

In the following result, we study the possibility of having the coexistence of the predator and prey populations.

Theorem 3.

Let ${\{(Z_n,{\tilde{Z}}_n)\}}_{n\in {\mathbb{N}}_0}$ be a predator–prey two-sex branching process. For any initial values $i,j>0$:

(i) 

If ${\tilde{\rho }}_1\tilde{m}\tilde{\alpha }<1$, then $P_{(i,j)}(Z_n\to \infty ,{\tilde{Z}}_n\to \infty )=0$.

(ii) 

If $min\{{\rho }_{2}m\alpha ,{\tilde{\rho }}_1\tilde{m}\tilde{\alpha }\}>1$, then $P_{(i,j)}(Z_n\to \infty ,{\tilde{Z}}_n\to \infty )>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 [14]), 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,45,46,47] or [48]). 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 [49,50] or [51]).

7. Discussion

In this paper a controlled two-sex branching process is introduced with the aim of modelling predator–prey interactions between populations with sexual reproduction for the first time in different situations. The generation-by-generation evolution of each species is also studied. We have focussed on the case where females and males form couples via a promiscuous sexual mating system in both species. 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 expected 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 [52]). 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 [53] or [54]). 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 [55]). 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(·)$ and $\tilde{s}(·)$ 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.