Research article Special Issues

Chance of extinction of populations in food chain model under demographic stochasticity

  • The extinction of different species from the earth is increasing at an alarming rate. So, assessment of probability of extinction of different important species in our ecosystem could help us to take proper conservation policy for those population whose chance of extinction is high. In this paper a method is developed to find the probability of extinction of populations in a general n-trophic food chain model under demographic stochasticity. The birth-death process is used to incorporate the demographic stochasticity and the necessary mathematical expressions are obtained. The theoretical finding is validated by numerical simulation for a two dimensional predator-prey system.

    Citation: Bapi Saha. Chance of extinction of populations in food chain model under demographic stochasticity[J]. Mathematical Biosciences and Engineering, 2019, 16(5): 3537-3560. doi: 10.3934/mbe.2019177

    Related Papers:

    [1] Aniello Buonocore, Luigia Caputo, Enrica Pirozzi, Amelia G. Nobile . A non-autonomous stochastic predator-prey model. Mathematical Biosciences and Engineering, 2014, 11(2): 167-188. doi: 10.3934/mbe.2014.11.167
    [2] Maoxiang Wang, Fenglan Hu, Meng Xu, Zhipeng Qiu . Keep, break and breakout in food chains with two and three species. Mathematical Biosciences and Engineering, 2021, 18(1): 817-836. doi: 10.3934/mbe.2021043
    [3] Claudio Arancibia–Ibarra, José Flores . Modelling and analysis of a modified May-Holling-Tanner predator-prey model with Allee effect in the prey and an alternative food source for the predator. Mathematical Biosciences and Engineering, 2020, 17(6): 8052-8073. doi: 10.3934/mbe.2020408
    [4] Xueqing He, Ming Liu, Xiaofeng Xu . Analysis of stochastic disease including predator-prey model with fear factor and Lévy jump. Mathematical Biosciences and Engineering, 2023, 20(2): 1750-1773. doi: 10.3934/mbe.2023080
    [5] Christian Cortés García, Jasmidt Vera Cuenca . Impact of alternative food on predator diet in a Leslie-Gower model with prey refuge and Holling Ⅱ functional response. Mathematical Biosciences and Engineering, 2023, 20(8): 13681-13703. doi: 10.3934/mbe.2023610
    [6] Lazarus Kalvein Beay, Agus Suryanto, Isnani Darti, Trisilowati . Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey. Mathematical Biosciences and Engineering, 2020, 17(4): 4080-4097. doi: 10.3934/mbe.2020226
    [7] Md Nazmul Hassan, Angela Peace . Mechanistically derived Toxicant-mediated predator-prey model under Stoichiometric constraints. Mathematical Biosciences and Engineering, 2020, 17(1): 349-365. doi: 10.3934/mbe.2020019
    [8] Dongmei Li, Tana Guo, Yajing Xu . The effects of impulsive toxicant input on a single-species population in a small polluted environment. Mathematical Biosciences and Engineering, 2019, 16(6): 8179-8194. doi: 10.3934/mbe.2019413
    [9] Jim M. Cushing . The evolutionarydynamics of a population model with a strong Allee effect. Mathematical Biosciences and Engineering, 2015, 12(4): 643-660. doi: 10.3934/mbe.2015.12.643
    [10] Adnan Sami, Amir Ali, Ramsha Shafqat, Nuttapol Pakkaranang, Mati ur Rahmamn . Analysis of food chain mathematical model under fractal fractional Caputo derivative. Mathematical Biosciences and Engineering, 2023, 20(2): 2094-2109. doi: 10.3934/mbe.2023097
  • The extinction of different species from the earth is increasing at an alarming rate. So, assessment of probability of extinction of different important species in our ecosystem could help us to take proper conservation policy for those population whose chance of extinction is high. In this paper a method is developed to find the probability of extinction of populations in a general n-trophic food chain model under demographic stochasticity. The birth-death process is used to incorporate the demographic stochasticity and the necessary mathematical expressions are obtained. The theoretical finding is validated by numerical simulation for a two dimensional predator-prey system.


    The rapid fall in the survivability of different important species is a major issue in ecology. The unplanned activity and indiscriminate use of biological resources lead to the unprecedented depletion in the abundance of species [1]. The detail study regarding conservation of different species is important since the loss of a strongly interactive species can lead to potential changes in ecosystem composition, structure, and diversity [2,3,4,5]. For example, the killing of large amount of great whales by industrial whaling makes changes in the krill-consumer dynamics in the Southern Ocean, and it is observed that, whaling is the main cause of megafaunal collapse in the North Pacific Ocean [6]. Moreover, the functional dynamics of coastal marine ecosystems worldwide have been extremely altered due to over fishing of large herbivores and predators [7]. It is frequently observed that the functional extinction of species interactions takes place well before the species themselves become extinct. In ocean ecosystem, large number of interactive species persists as rare adults, or small or juvenile individuals and do not interact as large adults. In addition, many strongly interactive species such as mammals are becoming extinct from the areas which they occupy one or two century ago [8]. The problem continues to persist since most conservation laws, such as the US Endangered Species Act (ESA) of 1973, fail to capture the effects of extinction of strongly interacting species. It is to be noted that the ESA does not include the inter-species interactions altogether in the criteria for recovery of endangered mammal species [4]. It only emphasizes the short-term, single-species, demographic viability in only a few areas. In reality, most of the current recovery plans do not call for increase in numbers of individuals or numbers of populations and geographic range [9,10]. Several works are also performed [11,12,13,14] to discuss this important issue. The persistence time of the population for a dynamical system is obtained in [15,16] which can be used to find the time to extinction of a population. However, as far as my knowledge goes the mathematical form to find the probability of extinction of a species belonging to any trophic level is still missing.

    Stochastic population models serve as fundamental mathematical framework in modern ecological theory and is applied in various research areas such as population ecology [17], epidemiology [18], conservation biology [19] etc. Among all these areas, epidemiology is a major platform where the knowledge of stochasticity is extensively used in [20,21,22] and [23]. It is observed that discrete and continuous time Markov chain models and stochastic differential equation models are used in many areas of population biology. The Kolmogorov differential equations (often termed as Master equation [24]) are used to compute the transition probability distributions when the rates of different possible transitions are provided (e.g. birth and death rates).

    It is frequently observed that to model a dynamical problem under stochastic environment, a system of Itostochastic differential equations is formed and studied. In many cases discrete stochastic models are developed by studying changes in the system components over a small interval of time such as birth-death process.

    The prey-predator models regulated by deterministic differential or difference equations are crucial in quantitative studies of the dynamics of natural populations. In general, the dynamical characteristics of the systems, such as stability of equilibrium points, stability of limit cycles etc. can be predicted uniquely in the deterministic set up. However, to apply these models in natural populations, it is more realistic to consider the stochastic dynamics of the system, where the randomness is generated either from variation among individual growth rates or due to environmental fluctuations. For such cases, the characteristics of the deterministic models are replaced by the probability statements. As for example, we may consider stochastic version of the logistic model, where the equilibrium point is replaced by the equilibrium distribution, which is represented by the cloud points around the equilibrium point rather than a fixed value [25]. The difference is that, in logistic growth, the carrying capacity is the only stable point, hence the population never reaches zero for positive value of intrinsic growth rate, whereas, in stochastic environment there is always a probability associated with the extinction of the population. Thus the stochastic model framework is more appealing in real life scenarios [26,27].

    The stability related issues of a predator-prey system is well discussed from the theoretical work of Lotka [28], Volterra [29] and the experimental work of Gause[30]. In many cases the dynamics involving only two interacting species may not capture the real world scenario. The behavior of a complex system may only be understood when the interactions among a large number of species are incorporated [31,32]. Thus, it should be noted that what amount of risk of extinction the interacting species face for a given initial population size of prey and predator, is an important aspect of population dynamics. Here the risk of extinction is quantified as the probability of reaching the extinction equilibrium, before reaching to a sufficiently large population size so that extinction can be avoided. Theoretically, this large population tends to [33,34]. Recently [34] used the stochastic differential equation model including demographic stochasticity and predicted the time to extinction and extinction probability for the Atlantic Herring populations. They calculated the mean time to extinction from the distribution of the sojourn time, the amount of time spent by the species at each population size given an initial abundance [35].

    In this paper, we explore the idea of obtaining the chance of extinction of species at any trophic level of the n trophic food chain model. Discrete stochastic model is used to study the changes in the system components over a small time interval to develop the probability of extinction. The theoretical findings are validated by numerical simulation.

    In next section (section 2), the background of the knowledge to be applied is discussed. In section 3, the main result of the paper is established with detail analysis followed by some case studies to understand the applicability of the findings. Section 4 is the numerical part followed by section 5 which deals with conclusion and future direction.

    In this section a two dimensional general predator prey system is considered as below,

    dxdt=b1(x)d1(x)g(x,y)dydt=αg(x,y)d2(y) (2.1)

    where x and y denote the sizes of prey and predator population at time t, respectively; b1(.) is the birth rate of the prey species; di(.)(i=1,2) are the natural death rates of the prey and predator populations respectively. The initial conditions are, x(0)0,y(0)0. The birth and death rates may be density dependent or density independent, described by the biological background of the two interacting species. The function g(x,y) denotes the inter-species interaction between the prey and predators, known as predator's functional response. α is the conversion rate of ingested prey into new predators and assumed to lie between 0 and 1.

    The stochastic formulation of the single species dynamics in absence of predators are studied by many authors [33,34]. In absence of predators, model (2.1) reduces to a single species model with birth and death rates b1(x) and d1(x) respectively. The logistic model under stochastic environment are studied in different ways in the literature [36,37,38,39,40]. In general, there are infinitely many choices of functions available for birth and death rates for different choices of parameters. The choices of the functions are mainly driven by the collected data on the species under investigation [41]. In addition, when we study real populations, generally some other information is also available such as, life expectancy. For example, the estimate of average life expectancy was used to obtain the estimates of intrinsic rate parameters[37] to model the dynamics of muskrat population in Netherlands and it help them to uniquely chose the value of the rate parameters.

    The generalized form of the model (2.1) in n-trophic food chain serves as the deterministic skeleton to formulate analogous stochastic models which account for the variation in births, deaths, transmission and recovery. In the stochastic formulation of the model (2.1) we assume X(t) and Y(t) are two discrete random variables representing the number of prey and predator populations respectively at time t and these two random variables take values in the state space {(x,y):x=0,1,M;y=0,1,N} where M and N are the maximum sustainable population sizes of prey and predator population respectively. Here we assume that the time step Δt is sufficiently small so that at most one change in the population size is possible in this small time gap, i. e. ΔX(t)=X(t+Δt)X(t){1,0,1} and ΔY(t)=Y(t+Δt)Y(t){1,0,1}. Thus the change in population sizes ΔX(t) and ΔY(t) preclude the possibility of more than one birth, death, or any other transformations in small time Δt. In fact, these events have probabilities of order (Δt)2. The infinitesimal probabilities associated with the birth and death of the populations for the model (2.1) is depicted in table 1. For example, (ΔX(t),ΔY(t))=(1,0) denotes the event of the birth of a prey and no change in predator abundance during time Δt. The probability of this event is obtained by multiplying the probability of one birth of prey (b1(x)Δt) with the probability of no death of predator ((1d2(y)Δt)). Similarly, the probabilities of the other events can be defined.

    Table 1.  Possible transitions in the predator-prey population sizes and the probabilities of their occurrence for the model (2.1).
    EventTransitionTransition probability
    Prey birth and no birth or death of predator(1,0)b1(x)Δt(1d2(y)Δt)
    One birth of predator due to predation on prey(1,1)αg(x,y)Δt
    No birth or death of prey and one death of predator(0,1)d2(y)Δt(1b1(x)Δt)
    Death of one prey due to intra-species competition or predation and no birth or death of predator(1,0)d1(x)Δt+(1α)g(x,y)Δt
    No birth or death of either population(0,0)1[b1(x)+αg(x,y)+d2(y)
    +d1(x)+(1α)g(x,y)]Δt

     | Show Table
    DownLoad: CSV

    We extend this idea for a tri-trophic food chain model given by,

    dxdt=b1(x)d1(x)g1(x,y)dydt=α1g1(x,y)g2(y,z)d2(y)dzdt=α2g2(y,z)d3(z) (2.2)

    where the transitions and the corresponding probabilities are given in the table 2.

    Table 2.  Possible transitions in the predator-prey population sizes and the probabilities of their occurrence for the tri-trophic food chain model (2.2).
    EventTransitionTransition probability (pi)
    Prey birth and no birth or death of predator(1,0,0)b1(x)Δt=p1
    One death of prey(1,0,0)((1α1)g1(x,y)+d1(x))Δt=p2
    One death of prey and one birth of 1st predator due to predation of 1st predator on prey(1,1,0)α1g1(x,y)Δt=p3
    No change of prey population, one death of 1st predator and one birth of 2nd predator due to predation of 2nd predator on the 1st predator(0,1,1)α2g2(y,z)Δt=p4
    No change of prey and 2nd predator and one death of 1st predator population(0,1,0)((1α2)g2(y,z)+d2(y))Δt=p5
    No change of prey and 1st predator and one natural death of 2nd predator(0,0,1)d3(z)Δt=p6
    No change in predator and prey population(0,0,0)1(p1+p2+p3+p4+p5+p6)

     | Show Table
    DownLoad: CSV

    We extend this further for an n trophic food chain model also as given below.

    The n trophic food chain model considered here is given by:

    dx1dt=b1(x1)g1(x1,x2)d1(x1)f1dx2dt=α1g1(x1,x2)g2(x2,x3)d2(x2)f2dx3dt=α2g2(x2,x3)g3(x3,x4)d3(x3)f3...dxndt=αn1gn1(xn1,xn)dn(xn)fn (2.3)

    Here the functions bi(xi)'s, di(xi)'s and gi(xi,xi+1) denote the birth rates, death rates and the interaction terms respectively.

    It is assumed that Xi(t)(i=1,2,,n) are the discrete random variables representing the number of prey and predator populations in the food chain. X(t) is the vector random variable which represents the population size at time t, that is X(t)=(X1(t),X2(t),,Xn(t)) and it takes values in the state space {x=(x1,x2,,xn):xi=0,1,Mi}. Mi is the maximum sustainable population size of the population belonging to ith trophic level. As before the time step Δt is taken sufficiently small so that only one change is possible in this small time interval, i. e. ΔXi(t)=Xi(t+Δt)Xi(t){1,0,1}. The infinitesimal probabilities of the birth death formulation of model (2.3) is depicted in table 3. For example, ΔX(t)=(1,0,,0) denotes the event of the birth of a prey and no change in the abundance of remaining species in time Δt. The probability of this event is (b1(x)Δt). Similarly probabilities of the other events are defined in table 3.

    Table 3.  Possible transitions in the predator-prey population sizes and the probabilities of their occurrence for the n-trophic food chain model (2.3).
    EventTransitionTransition probability (pi)
    One birth of prey and no change in population size of the remaining species(1,0,,0)b1(x)Δt=p1
    One birth of 1st predator due to predation on the prey population(1,1,0,,0)α1g1(x1,x2)Δt=p2
    One death of prey population and no change in the size of the populations of other species(1,0,0,,0)d1(x)Δt+(1α1)g1(x1,x2)Δt=p3
    One birth of 2nd predator due to predation on 1st predator(0,1,1,0,,0)α2g2(x2,x3)Δt=p4
    One death of 1st predator and no change in the population size in other species(0,1,0,,0)d2(x2)Δt+(1α2)g2(x2,x3)Δt=p5
    One birth of (n1)th predator due to predation on (n2)th predator(0,,0,1,1)αn1gn1(xn1,xn)Δt=p2n2
    One death of (n2)th predator and no change in the population size of other species(0,,0,1,0)dn1(xn1)Δt+(1αn1)gn1(xn1,xn)Δt=p2n2
    One death of (n1)th predator and no change in the population size of other species(0,,0,1)dn(xn)Δt=p2n
    No change in the population size in any of the species in the food chain(0,0,,0)1Σ2ni=1pi

     | Show Table
    DownLoad: CSV

    It should be noted that the total number of transitions for an n trophic food chain model is 2n+1. The expected values of the changes in small time interval Δt may be obtained as follows.

    E(ΔX1)=1.b1(x1)Δt+(1)(α1g1(x1,x2)+d1(x1)+(1α1)g1(x1,x2))Δt=(b1(x1)d1(x1)g1(x1,x2))Δt[see table 3]=dx1dtΔt (2.4)

    Similarly, we can find

    E(ΔXi)=1.αi1gi1(xi1,xi)+(1)(αigi(xi,xi+1)+(1αi)gi(xi,xi+1)+di(xi))Δt=(αi1gi1(xi1,xi)gi(xi,xi+1)di(xi))Δtfori=2,3,,n1=dxidtΔt (2.5)

    Finally we have,

    E(ΔXn)=(αn1gn1(xn1,xn)dn(xn))Δt=dxndtΔt (2.6)

    Now we will find the expressions for variance and co-variance terms for the changes.

    We have,

    E(ΔX2i)=12αi1gi1(xi1,xi)+(1)2(gi(xi,xi+1)+di(xi))Δt=(αi1gi1(xi1,xi)+gi(xi,xi+1)+di(xi))Δt (2.7)
    fori=2,3,,n1 (2.8)
    E(ΔX21)=(b1(x1)+d1(x1)+g1(x1,x2))Δt (2.9)
    E(ΔX2n)=(αn1gn1(xn1,xn)+dn(xn))Δt (2.10)

    Now it can be shown that,

    E(ΔxiΔxi1)=αi1gi1(xi1,xi)Δtfori=2,3,,n1 (2.11)

    It should be noted that, here the randomness in population growth rate is assumed to be demographic stochasticity which is the chance variation in the number of individual births and deaths and usually modeled by birth and death process. So the above model does not depict the random variations in the environmental conditions although the environmental fluctuation has significant effect on the population dynamics, for example it may have effect on the survival and reproduction rates [42].

    In this section the expression for the probability of extinction of a particular species population in the given food chain model is derived. Let Tx(s) be the first passage time of the random variable X to the state s. This is also known as persistent time or the first exit time in engineering literature. Thus the event Tx(0)<Tx(b) represents the event that starting with population size x, the species goes to extinction before reaching to the population size b.

    An n dimensional food chain model as given below is considered here.

    dx1dt=b1(x1)g1(x1,x2)d1(x1)f1dx2dt=α1g1(x1,x2)g2(x2,x3)d2(x2)f2dx3dt=α2g2(x2,x3)g3(x3,x4)d3(x3)f3...dxndt=αn1gn1(xn1,xn)dn(xn)fn (3.1)

    Suppose we are interested to find the chance of extinction of ith species. We start with the 1st species (here the prey population) and subsequently do the same for ith species.

    Define,

    u1(x1,x2,,xn)=P[Tx1(0)<Tx1(M1) (3.2)
    |X1(0)=x1,X2(0)=x2,,Xn(0)=xn] (3.3)

    In this case it should be noted that, according to the definition u1 is not a random variable. We now derive the pde for u1. It can be shown that,

    u1(x1,x2,,xn)=P[Tx1(0)<Tx1(M1)|X1(0)=x1,,Xn(0)=xn]=E[u1(x1+Δx1,x2+Δx2,,xn+Δxn)]=E[u1(x1,x2,,xn)+Δx1u1x1++Δxnu1xn+ijΔxiΔxj2u1xixj][ignoring higher order terms]0=E[Δx1]u1x1+E[Δx2]u1x2++E[Δxn]u1xn+ (3.4)
    E[(Δx1)2]2u1x21++E[(Δxn)2]2u1x2n+2E[Δx1Δx2]2u1x1x2++2E[Δxn1Δxn]2u1xn1xn0=u1x1dx1dt+u1x2dx2dt++u1xndxndt+(b1(x1)+d1(x1)+g1(x1,x2))2u1x21+(αn1gn1(xn1,xn)+dn(xn))2u1x2nα1g1(x1,x2)2u1x1x2αn1gn1(xn1,xn)2u1xn1xn (3.5)

    see (2.4), (2.5), (2.6), (2.7), (2.11)

    Now using the expressions of E(ΔXi), E(ΔX2i) and E(ΔXiΔXi1) and noting that

    u1t=u1x1dx1dt+u1x2dx2dt++u1xndxndt (3.6)

    we get,

    u1t=α1g1(x1,x2)2u1x1x2+α2g2(x2,x3)2u1x2x3++αn1gn1(xn1,xn)2u1xn1xn(b1(x1)+g1(x1,x2)+d1(x1))2u1x21(α1g1(x1,x2)+g2(x2,x3)+d2(x2))2u1x22(αn1gn1(xn1,xn)+dn(xn))2u1x2n (3.7)

    Theoretically, probability of extinction of the prey species can be obtained by letting M1. When the system is in steady state u1t=0.

    So, in steady state the above pde takes the form,

    0=α1g1(x1,x2)2u1x1x2+α2g2(x2,x3)2u1x2x3++αn1gn1(xn1,xn)2u1xn1xn(b1(x1)+g1(x1,x2)+d1(x1))2u1x21(α1g1(x1,x2)+g2(x2,x3)+d2(x2))2u1x22(αn1gn1(xn1,xn)+dn(xn))2u1x2n (3.8)

    The boundary conditions are given by,

    u1(0,x2,,xn)=1 (3.9)
    u1(M1,x2,,xn)=0 (3.10)
    u1(x1,0,,xn)=f1(x1) (3.11)
    u1(x1,,Mi,xn)=ϕi(x1,,xi1,xi+1,,xn)say (3.12)

    If the system starts with 0 population for 1st predator, the whole system becomes a single species dynamical system consisting of the prey population only. Because in this case all the predator population will die with probability 1. Hence, the said system will be governed by single species model and there will be no interaction term. In that case u1(x,,0)=f1(x). The expression for f1 is given by [35],

    f1(x)=M1xexp[ϕ(u)]duM10exp[ϕ(u)]duwhereϕ(u)=2b1(u)d1(u)b1(u)+d1(u)du

    It can be shown that ϕi satisfies (3.4) by setting Δx2i=0 and ΔxkΔxi=0 (k=1,2,,n,ki) because we assume that when the ith species reaches the state Mi, the small change in the population size at ith trophic in time Δt will be negligible compared to Mi.

    To find the chance of extinction of the ith species we define,

    ui(x1,x2,,xn)=P[Txi(0)<Txi(Mi)|X(0)=x,Txi1(0)>Txi1(Mi1)]ifxi,xi10=1ifxi=0orxi1=0 (3.13)

    where X=(X1,X2,,Xn) and x=(x1,x2,,xn).

    It should be noted that, in order to find the chance of extinction for ith(i=2,3,,n)population, the condition that the (i1)th species do not die out (Txi1(0)>Txi1(Mi1)) is imposed. This condition is necessary because the loss of a species leads to the extinction of other dependant species [1]. ui(x) represents probability of extinction of ith species when Mi. Now

    ui(x1,x2,,xn)=P[Txi(0)<Txi(Mi)|X(0)=x,Txi1(0)>Txi1(Mi1)]=P[Txi(0)<Txi(Mi),Txi1(0)>Txi1(Mi1)|X(0)=x]P[Txi1(0)>Txi1(Mi1)|X(0)=x]ui(x)=wi1(x)1ui1(x)ifxi,xi10i = 2, 3, ,n=1ifxiorxi1=0 (3.14)

    where,

    wi1(x)=P[Txi(0)<Txi(Mi),Txi1(0)>Txi1(Mi1)|X(0)=x] (3.15)

    Following the same idea as applied in case of deriving the pde for u1, it can be shown that wi1 satisfies the same pde (3.7) as u1 does with the following boundary conditions.

    wi1(0,x2,,xn)=0 (3.16)
    wi1(x1,x2,,Mi,,xn)=0 (3.17)
    wi1(x1,x2,,Mk,,xn)=ψki1(x1,,xk1,xk+1,,xn) (3.18)
    forki (3.19)

    This ψki1 satisfies the pde obtained by setting Δx2k=0 and ΔxkΔxj=0 in (3.4) for j=2,,n and jk.

    So far the discussion revolves around the general n-trophic food chain model. In this section two particular cases viz. two dimensional and three dimensional food chain model for better understanding are discussed.

    First, a two dimensional food chain model (2.1) is considered.

    The probability of extinction for theses types of model is also provided in [43].

    We define:

    u(x,y)=P{Tx(0)<Tx(M)|X(0)=x,Y(0)=y} (3.20)

    where Tx(.) is defined earlier. u(x,y) represents the probability of extinction when the dynamical system start from the state (x,y) allowing M [33]. Now,

    u(x,y)=ΣΔxΣΔyP{Tx(0)<Tx(M),X(0+Δt)=x+Δx,Y(0+Δt)=y+Δy|X(0)=x,Y(0)=y} (3.21)

    From [35], it can be shown that

    u(x,y)=E[u(X(Δt),Y(Δt))|X(0)=x,Y(0)=y] (3.22)

    Now,

    u(x+Δx,y+Δy)=u(x,y)+Δxux+Δyuy+12(Δx)22ux2+ΔxΔy2uxy+122uy2(Δy)2

    Taking expectation on both sides and using (3.22) we get,

    u(x,y)=u(x,y)+E(Δx)ux+E(Δy)uy+12E(Δx)22ux2+E(ΔxΔy)2uxy+122uy2E(Δy)2 (3.23)

    Using the equation (2.1), we obtain the following equation

    0=[b1(x)g(x,y)d1(x)]ux+[αg(x,y)d2(y)]uy+12[b1(x)+g(x,y)+d1(x)]2ux2αg(x,y)2uxy+12[αg(x,y)+d2(y)]2uy2

    Now replacing the coefficients of ux and uy by dxdt and dydt respectively, by virtue of the equation (2.1), and noting that, ut=uxdxdt+uydydt, we obtain the following partial differential equation

    0=ut+12[b1(x)+g(x,y)+d1(x)]2ux2αg(x,y)2uxy+12[αg(x,y)+d2(y)]2uy2

    Therefore the above expression can be written as

    ut=12[b1(x)+g(x,y)+d1(x)]2ux2+αg(x,y)2uxy12[αg(x,y)+d2(y)]2uy2 (3.24)

    In steady state the pde (3.24) takes the form:

    12[b1(x)+g(x,y)+d1(x)]2ux2αg(x,y)2uxy+12[αg(x,y)+d2(y)]2uy2=0 (3.25)

    Lemma 3.1. The pde (3.25) is elliptic in nature.

    Proof. The standard form of a second order partial differential equation is

    A2ux2+B2uxy+C2uy2+Dux+Euy+Fu+G=0 (3.26)

    In this case A=12[b1(x)+g(x,y)+d1(x)], B=αg(x,y) and C=12[αg(x,y)+d2(y)]. Therefore, B24AC=α(1α)g2(x,y)d2(y)(b1(x)+d1(x))g(x,y)(αb1(x)+αd1(x)+d2(y))<0 since g(x,y),b1(x),d1(x),d2(y) are all positive and α<1.

    Hence the pde (3.25) is elliptic.

    To compute the first passage probabilities given the size of prey and predator, the above equation is solved numerically. We now proceed to develop appropriate boundary conditions for the above differential equation. We have,

    u(x,y)=P{Tx(0)<Tx(M)|X(0)=x,Y(0)=y}

    Therefore

    u(0,y)=P{Tx(0)<Tx(M)|X(0)=0,Y(0)=y}=1 (3.27)

    because, if the initial prey population is 0, the event Tx(0)<Tx(M) is a certain event.

    Also,

    u(x,0)=P{Tx(0)<Tx(M)|X(0)=x,Y(0)=0}=P{Tx(0)<Tx(M)|X(0)=x}=f(x)say (3.28)

    The function f(x) is given by ([35])

    f(x)=Mxexp[ϕ(u)]duM0exp[ϕ(u)]duwhere (3.29)
    ϕ(u)=2b(u)d(u)b(u)+d(u)du (3.30)

    Now

    u(M,y)=P{Tx(0)<Tx(M)|X(0)=M,Y(0)=y}=0,becauseP{Tx(0)<Tx(M)|X(0)=M}=0 (3.31)
    u(x,N)=P{Tx(0)<Tx(M)|X(0)=x,Y(0)=N}=g(x)(say) (3.32)

    The differential equation for g(x)=u(x,N) can be derived by setting the coefficients of 2uy2 and 2uxy in (3.23) to be 0 as N is assumed to be the maximum possible size of predator population. In addition to this ut=0 for stationary state. Thus,

    d2gdx2=0 (3.33)

    The boundary conditions are g(0)=u(0,N)=1,g(M)=u(M,N)=0. Using this, the solution of the above differential equation is,

    g(x)=u(x,N)=1xM (3.34)

    Remark 1. We have,

    limMu(x,N)=limM(1xM)=1 (3.35)

    This indicates that if the predator population is at its maximum level, the extinction of the prey population is inevitable if the prey population is sufficiently small compared to the maximum size (M) of the prey population.

    Now we will develop the method of finding the probability of extinction of predator population. We define

    v(x,y)=P{Ty(0)<Ty(N)|Tx(0)>Tx(M),X(0)=x,Y(0)=y}ifx,y0=1ifx=0ory=0 (3.36)

    The condition Tx(0)>Tx(M) is imposed since otherwise the extinction of predator species is certain. We can write,

    v(x,y)=P{Ty(0)<Ty(N),Tx(0)>Tx(M)|X(0)=x,Y(0)=y}P(Tx(0)>Tx(M))=w(x,y)1u(x,y) (3.37)

    where

    w(x,y)=P{Ty(0)<Ty(N),Tx(0)>Tx(M)|X(0)=x,Y(0)=y} (3.38)

    It can be shown that (in a similar way as in the case of u) w satisfies the same pde as u with the boundary conditions given next.

    We have,

    w(0,y)=v(0,y)(1u(0,y))=0(u(0,y)=1) (3.39)

    w(M,y) will be obtained as in case of u(x,N) and is found to be,

    w(M,y)=1yN (3.40)

    Remark 2. We have

    v(M,y)=w(M,y)1u(M,y)=w(M,y)10=w(M,y) (3.41)

    Now it should be noted that,

    limNw(M,y)=limN(1yN)=1 (3.42)

    The implication of this result is that if the prey population size is at its maximum, the predator population will reach state 0 before it reaches its maximum (i. e. N) if N is very large compared to the initial population size of predator.

    The remaining boundary conditions are,

    w(x,0)=v(x,0)(1u(x,0))=1u(x,0)sincev(x,0)=1 (3.43)
    w(x,N)=0 (3.44)

    In many cases tri-trophic food chain model is appropriate to capture more information for a complex system [44]. In order to find the chance of extinction of a population at any trophic level, we can extend the above idea for tri-trophic food chain model also. Here a tri-trophic food chain model of the form given below is considered,

    dxdt=b1(x)g1(x,y)d1(x)dydt=α1g1(x,y)g2(y,z)d2(y)dzdt=α2g2(y,z)d3(z) (3.45)

    Similar procedure may be used to define,

    u1(x,y,z)=P{Tx(0)<Tx(M)|X(0)=x,Y(0)=y,Z(0)=z} (3.46)

    We can follow the same argument as in the case of general n-trophic food chain model to show that,

    u1t=12[b1(x)+g1(x,y)+d1(x)]2u1x212[α1g1(x,y)+g2(y,z)+d2(y)]2u1y212[α2g2(y,z)+d3(z)]2u1z2+α1g1(x,y)2u1xy+α2g2(y,z)2u1yz (3.47)

    The six boundary conditions on the six surfaces are given by,

    u1(0,y,z)=1u1(M1,y,z)=0u1(x,0,z)=f(x)u1(x,M2,z)=f1(x,z)sayu1(x,y,0)=f2(x,y)sayu1(x,y,M3)=f3(x,y)say (3.48)

    The function f1(x,z) can be found in the following way:

    Since M2 is assumed to be the maximum possible population size of 1st predator, the decrease of it by one individual is negligible compared to M2. Hence we may set Δx2=0 and (Δx2)2=0 in the general equation (3.4) for 1st predator.

    This will give the following differential equation in stationary state,

    [b1(x)+g1(x,y)+d1(x)]2f1x2+[α2g2(y,z)+d3(z)]2f1z2=0 (3.49)

    f1 satisfies the boundary conditions,

    f1(M1,z)=u1(M1,M2,z)=0f1(x,M3)=1xM1f1(0,z)=1f1(x,0)=u1(x,M2,0)=1xM1 (3.50)

    In case of u1(x,y,0), the system is free from 3rd predator and hence it is equivalent to 2dimensional prey-predator system. This suggest that u1(x,y,0)=u(x,y)=f2(x,y).

    Following the argument in order to derive the differential equation for f3, it can be shown that,

    (b1(x)+g1(x,y)+d1(x))2f3x2+(α1g1(x,y)+g2(y,M3)+d2(y))2f3y22α1g1(x,y)2f3xy=0 (3.51)

    The boundary conditions for (3.51) are,

    f3(M1,y)=0f3(0,y)=1f3(x,M2)=1xM1f3(x,0)=f(x)(as given in the previous section) (3.52)

    In case of 1st predator population we define,

    u2(x,y,z)=P{Ty(0)<Ty(M2)|Tx(0)>Tx(M1),X(0)=x,Y(0)=y,Z(0)=z}ifx,y0=1ifx=0ory=0 (3.53)

    As in the case of two dimensional model, the expression for u2 can be written in the form,

    u2(x,y,z)=w1(x,y,z)1u1(x,y,z)ifx,y0=1ifx=0ory=0 (3.54)

    where,

    w1(x,y,z)=P[Ty(0)<Ty(M2),Tx(0)>Tx(M1)|X(0)=x,Y(0)=y,Z(0)=z] (3.55)

    w1 satisfies the same pde as u1 with suitable boundary conditions given below.

    w1(0,y,z)=0. From the definition of w1 we have,

    w1(x,0,z)=P[Ty(0)<Ty(M2),Tx(0)>Tx(M1)|X(0)=x,Y(0)=0,Z(0)=z]=P[Tx(0)>Tx(M1)|X(0)=x,Y(0)=0,Z(0)=z]×P[Ty(0)<Ty(M2)|Tx(0)>Tx(M1),X(0)=x,Y(0)=0,Z(0)=z]=(1u1(x,0,z))u2(x,0,z)=(1u1(x,0,z)).1 (3.56)

    Hence we have w1(x,0,z)=1u1(x,0,z)=1f(x). It can be shown that w1(x,M2,z)=0. In order to find w1(x,y,0) we follow the following procedure.

    Let w1(x,y,0)=ϕ(x,y)(say). It can be shown that ϕ(x,y) satisfies the pde (3.3.2) with the coefficients of 2ϕz2, 2ϕyz to be 0. Hence ϕ(x,y) satisfies at steady state,

    0=[b1(x)+α1g(x,y)+d1(x)]2ϕx2+[α1g1(x,y)+g2(y,z)+d2(y)]2ϕy22α1g1(x,y)2ϕxy (3.57)

    ϕ(x,y) satisfies the following boundary conditions. ϕ(x,0)=w1(x,0,0)=1f(x) using the expression for w1(x,0,z) given above. ϕ(x,M2)=0, ϕ(0,y)=0 and ϕ(M1,y)=1yM2.

    Similarly, w1(x,y,M3)=ψ(x,y) can be shown to satisfy

    0=[b1(x)+α1g1(x,y)+d1(x)]2ψx2+[α1g1(x,y)+g2(y,M3)+d2(y)]2ψy22α1g1(x,y)2ψxy (3.58)

    The boundary conditions for ψ are:

    ψ(x,0)=1f(x)ψ(x,M2)=0ψ(0,y)=u2(0,y,M3)(1u1(0,y,M3))=u2(0,y,M3)(11)=0ψ(M1,y)=0w1(0,y,z)=0w1(M1,y,z)=ξ(y,z)(say) (3.59)

    which can be obtained in similar way.

    In continuation of this process of formulating the probability of extinction we can find the same (u3) for the 2nd predator also.

    In this section the theoretical findings are validated by using the following simple two dimensional Lotka-Volterra food chain model.

    dxdt=rx(1xK)xydydt=αxydy (4.1)

    Here r is the intrinsic growth rate of the prey population, K is the carrying capacity, α is the conversion coefficient and d is the natural death rate of the predator population. If we compare the model (4.1) with the model (2.1), we may consider the following birth-death and interaction terms.

    b1(x)=(r+1)xr2Kx2 (4.2)
    d1(x)=x+r2Kx2 (4.3)
    g(x,y)=xy (4.4)
    d2(y)=dy (4.5)

    The probability of extinction of prey population (u) and predator population (v) are generated numerically taking help of finite difference method using the procedure as given in section (3.3.1). Sufficiently large values of M and N are considered to have better approximation for probability of extinction. The surface of u(x,y) and v(x,y) are given in the figures 1 and 2.

    Figure 1.  The plot of u vs. the initial population sizes of prey and predator for the model (4.1). The parameters are r=0.2, K=20, α=0.5, d=0.1.
    Figure 2.  The plot of v vs. the initial population sizes of prey and predator for the model (4.1). The parameters are r=0.2, K=20, α=0.5, d=0.1.

    One can notice in the figure 1 that the probability of extinction is high in the region where the initial population sizes of both prey and predator population are low and vice-versa. Note that it is natural for any species to be prone to extinction if its initial population size is very small. This natural phenomena is captured in this proposed method. The solution set of u(x,y) and v(x,y) provide the probability of extinction of prey and predator population respectively.

    To extend the applicability of this method or to understand the theoretical findings in a better way, we also consider the effect of initial population size of either prey or predator separately on the extinction probability. We observe that if the prey population size increases, the probability of extinction of prey population diminishes for any fixed initial population size of predator(Figure 3). The interesting result we observe that, if the initial size of predator increases, the chance of extinction of prey decreases until the initial size of predator population reaches a certain value. The further increase in the initial size of predator increases the chance of extinction of prey population (see Figure 3 and 4). Thus the presence of a moderate amount of predator in the system is important to avoid extinction of prey which in turn helps the predator to survive also.

    Figure 3.  The plot of u vs. the initial population sizes of prey for three different initial population size of predator. The parameters are r=0.2, K=20, α=0.5, d=0.1.
    Figure 4.  The plot of u vs. the initial population sizes of predator for three different initial population size of prey. The parameters are r=0.2, K=20, α=0.5, d=0.1.

    In figure 5, we observe that if the initial prey population increases, the chance of extinction (v here) of predator decreases rapidly up to a certain value of initial prey. The change in the extinction probability of predator is not significant for further increase in the prey population and when initial prey population size exceeds a certain value, the chance of extinction of predator increases. Here also we observe that if a food chain model start with very high abundance of initial prey population predator population may become extinct with high probability. We also observe that for any fixed value of initial prey population, the higher the initial predator population, lower is the chance of extinction of it.

    Figure 5.  The plot of v vs. the initial population sizes of prey for three different initial population size of predator. The parameters are r=0.2, K=20, α=0.5, d=0.1.

    Indiscriminate use of biological resources along with unplanned activity of human [45] has adverse effect on ecosystem. As a result, a good number of rare and endangered species are becoming extinct very rapidly. If this continues to happen it can be predicted that the species in the ecosystem will be extinct in an unprecedented rate [46,47]. In this context, the assessment of extinction threats of different rare and endangered species are needed so that corresponding safeguard can be provided.

    In this paper we have formulated the technique to find the probability of extinction for an n-trophic food chain model. This process can be used to find the extinction probability of species belonging to any trophic level of the food chain and proper conservation policy can be implemented in advance. This strategy may help significantly to account for the extinction threat of species in widely recognized bio-diverse places like Atlantic Forest [48], Sumatra and Madagascar etc. [49]. It is observed that, the variables determining the extinction probabilities follow few partial differential equations some of which are very complex in nature and can not be solved explicitly. The complexity increases with the increase in the dimension of food chain. So numerical method may be used to find the solution of those differential equations to get the probability of extinction.

    In this paper we consider only demographic variability. The impact of environmental noise on any food chain is also significant. It will be interesting and more appropriate if this work is extended incorporating environmental noise also.

    I am thankful to Govt. College of Engg. & Textile Technology, Berhampore for this research. I also express my thank to fellow researchers with whom I discussed and take suggestions for this work and to the learned reviewers whose suggestions help me immensely to improve the manuscript.

    The authors declared they have no competing interests.



    [1] B. Ebenman, R. Law and C. Borrvall, Community viability analysis: The response of ecological
    [2] M. E. Soulé and J. Terborgh, Protecting nature at regional and continental scales: A conservation biology program for the new millennium, BioScience, 49 (1999), 809–817.
    [3] O. J. Schmitz, P. A. Hamback and A. P. Beckerman, Trophic cascades in terrestrial systems: A review of the effects of carnivore removal on plants, Am. Nat., 155 (2000), 141–153.
    [4] M. E. Soulé, J. Estes, J. Berger, et al., Ecological effectiveness: Conservation goals for interactive species, Conserv. Biol., 17 (2003), 1238–1250.
    [5] L. Oksanen and T. Oksanen, The logic and realism of the hypothesis of exploitation ecosystems, Am. Nat., 118 (2000), 240–261.
    [6] A. M. Springer, J. E. Estes, G. B. van Vliet, et al., Sequential megafaunal collapse in the North Pacific Ocean: An ongoing legacy of industrial whaling?, Proceed. Nat. Aca. Sci., 100 (2003), 12223–12228.
    [7] J. B. C. Jackson, M. X. Kirby, W. H. Berger, et al., Historical over fishing and the recent collapse of coastal ecosystems, Science, 293 (2001), 629–638.
    [8] A. S. Laliberte and W. J. Ripple, Range contractions of North American carnivores and ungulates, BioScience, 54 (2004), 123–138.
    [9] T. H. Tear, J. M. Scott, P. H. Hayward, et al., Recovery plans and the endangered species act:Are criticisms supported by data?, Conserv. Biol., 9 (1995), 182–195.
    [10] D. Jennings, South Florida multi-species recovery plan. vero beach (FL): US fish and wildlife service, (1999).
    [11] I. Nijs and I. Impens, Biological diversity and probability of local extinction of ecosystems, Funct. Ecol., 14 (2000), 46–54.
    [12] C. Borrvall, B. Ebenman and T. Jonsson, Biodiversity lessens the risk of cascading extinction in model food webs, Ecol. Lett., 3 (2000), 131–136.
    [13] Y. Ma and Q. Zhang, Stationary distribution and extinction of a three-species food chain stochastic model, Transact. A. Razmadze Math. Inst., 172 (2018), 251–264.
    [14] J. Baumsteiger and P. B. Moyle, Assessing Extinction, BioScience, 67 (2017), 357–366.
    [15] E. J. Allen, Stochastic differential equations and persistence time for two interacting populations, Dynam. Cont. Discret. Impul. Systems, 5 (1999), 271–281.
    [16] E. J. Allen, L. J. S. Allen and H. Schurz, A comparison of persistence-time estimation for discrete and continuous population models that include demographic and environmental variability, Math. Biosci., 196 (2005), 14–38.
    [17] M. J. Keeling, Multiplicative moments and measures of persistence in ecology, J. Theor. Biol., 205 (2000), 269–281.
    [18] I. Krishnarajah, A. Cook, G. Marion, et al., Novel moment closure approximations in stochastic epidemics, Bullet. Math. Biol., 67 (2005), 855–873.
    [19] R. Lande, S. Engen and B-E. Sæther, Stochastic Population Dynamics in Ecology and Conserva- tion, Oxford University Press, Oxford, 2003.
    [20] Y. Cai, Y. Kang, M. Banerjee, et al., A stochastic SIRS epidemic model with infectious force under intervention strategies, J. Differ. Equat., 259 (2015), 7463–7502.
    [21] Y. Cai, J. Jiao, Z. Gui, et al., Environmental variability in a stochastic epidemic model, Appl. Math. Comput., 329 (2018), 210–226.
    [22] W. Guo, Y. Cai, Q. Zhang, et al., Stochastic persistence and stationary distribution in an SIS epidemic model with media coverage, Physica A, 492 (2018), 2220–2236.
    [23] B. Yang, Y. Cai, K. Wang, et al., Global threshold dynamics of a stochastic epidemic model incorporating media coverage, Adv. Differ. Equat., 462 (2018).
    [24] N. T. J. Bailey, Element of Stochastic Processes with applications to the natural sciences, Wiley, New York, 1964.
    [25] R. M. May, Stability in randomly fluctuating versus deterministic environments, Am. Nat., 107 (1973), 621–650.
    [26] E. Renshaw, Modelling Biological Populations in Space and Time, Cambridge University Press, Cambridge, 1993.
    [27] E. Renshaw, Stochastic population processes: analysis, approximations, simulations, Oxord University Press, New York, 2011.
    [28] A. J. Lotka, Elements of physical biology, Science Progress in the Twentieth Century (1919-1933), 21 (1926), 341–343.
    [29] V. Volterra, Lecons sur la thorie mathmatique de la lutte pour la vie, Gauthier-Villars, Paris, 1931.
    [30] G. F. Gause, The Struggle for Existence, 1934.
    [31] S. L. Pimm, Food webs, Chapman and Hall, New York, USA, 1982.
    [32] R. T. Paine, Food webs: road maps of interaction or grist for theoretical development?, Ecology, 69 (1988), 1648–1654.
    [33] B. Dennis, Allee effects in stochastic populations, Oikos, 96 (2002), 389–401.
    [34] B. Saha, A. R. Bhowmick, J. Chattopadhyay, et al., On the evidence of an allee effect in herring populations and consequences for population survival: A model-based study, Ecol. Modell., 250 (2013), 72–80.
    [35] S. Karlin and H. M. Taylor, A second course in stochastic process, Oxord University Press, New York, 1981.
    [36] J. H. Matis and T. R. Kiffe, On approximating the moments of the equilibrium distribution of a stochastic logistic model, Biometrics, 52 (1996), 980–991.
    [37] J. H. Matis, T. R. Kiffe and P. R. Parthasarathy, On the cumulants of population size for the stochastic power law logistic model, Theor. Popul. Biol., 53 (1998), 16–29.
    [38] J. H. Matis and T. R. Kiffe, Effects of immigration on some stochastic logistic models: A cumulant truncation analysis, Theor. Popul. Biol., 56 (1999), 139–161.
    [39] I. Nåsell, Extinction and quasi-stationarity in the Verhulst logistic model, J. Theor. Biol., 211 (2001), 11–27.
    [40] I. Nåsell, Moment closure and the stochastic logistic model, Theor. Popul. Biol., 63 (2003), 159–168.
    [41] A. R. Bhowmick, B. Saha, S. Ray, J. Chattopadhyay and S. Bhattacharya, Cooperation in species: Interplay of population regulation and extinction through global population dynamics database, Ecol. modell., 312 (2015), 150–165.
    [42] R. Lande, Risks of population extinction from demographic and environmental stochasticity and random catastrophes, Am. Nat., 142 (2004), 911–927.
    [43] B. Saha, On extended growth model, functional response and related mathematical issues: Il- lustration through experimental, history and simulated data, Ph.D thesis, University of Calcutta, 2015.
    [44] A. Hastings and T. Powell, Chaos in a three species food chain, Am. Nat., 72 (1991), 896–903.
    [45] J. W. Bull and M. Maron, How humans drive speciation as well as extinction, Proc. R. Soc. B, 283 (2016).
    [46] S. L. Pimm, J. L. Gittleman and T. Brooks, The future of biodiversity, Science, 269 (1995), 347–350.
    [47] O. E. Sala, Global biodiversity scenarios for the year 2100, Science, 287 (2000), 1770–1774.
    [48] T. Brooks, J. Tobias and A. Balmford, Deforestation and bird extinctions in the Atlantic forest, Anim. Conserv., 2 (1999), 211–222.
    [49] N. Ocampo-Pe˜ nuela, C. N. Jenkins, V. Vijay, et al., Incorporating explicit geospatial data shows more species at risk of extinction than the current Red List, Sci. Adv., 2 (2016).
  • This article has been cited by:

    1. Bapi Saha, Rupak Bhattacharjee, Debasish Majumder, 2021, Chapter 5, 978-981-15-6886-2, 43, 10.1007/978-981-15-6887-9_5
    2. Sourav Kumar Sasmal, Yasuhiro Takeuchi, Editorial: Mathematical Modeling to Solve the Problems in Life Sciences, 2020, 17, 1551-0018, 2967, 10.3934/mbe.2020167
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4855) PDF downloads(638) Cited by(2)

Figures and Tables

Figures(5)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog