# Birth Rate Effects on an Age-Structured Predator-Prey Model with Cannibalism in the Prey

Authors: Francisco J. Solis; Roberto A. Ku-Carrillo
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

We develop a family of predator-prey models with age
structure and cannibalism in the prey population. It consists of systems
ofm ordinary differential equations, where m is a parameter associated with
new proposed prey birth rates. We discuss how these new birth rates give
the required flexibility to produce differential systems with well-behaved
solutions. The main feature required in these models is the coexistence
among the involved species, which translates mathematically into stable
equilibria and periodic solutions. The search for such characteristics is
based on heuristic predation functions that account for cannibalism in
the prey.


---

## Body

## 1. Introduction

The modeling of predator-prey populations is an intensive research area where a wide variety of models are proposed to describe and to analyze different mechanisms such as control of populations, oscillations of populations, delay maturation time effects, and pattern formation in predator-prey models with diffusion: see, for example, [1–10]. One relevant feature in these models is the predator-prey coexistence, which is of vital importance to obtain realistic models of predator-prey interactions. Mathematically, a first type of coexistence is established by stable equilibria which are the most basic solutions to the models. A second type of coexistence is given by stable periodic solutions and this is the most common accepted type in predator-prey ecological systems; however, it is technically more difficult to achieve.Our first goal is to propose and analyze models derived from the classical approach suggested by Gurtin and MacCamy [5, 6] for predator-prey population models with age structure. Such approach has a starting point given by the McKendrick-Von Foester equation [11] and the generated models involve integrodifferential equations. By introducing appropriate birth rates, it is possible to obtain models given by systems of ordinary differential equations. An important characteristic of our proposed models is the consideration of the cannibalism effect on the prey, since several studies have reported the stabilizing consequences of cannibalism in predator-prey systems [12, 13]. There are numerous examples in nature where cannibalism is present, for example, the confused flour beetle (Tribolium confusum) [14], the red flour beetle (Tribolium castaneum) [15], and the Dungeness crab (Metacarcinus magister) [16].Our second goal is to analyze the relationship between the existence of periodic solutions and the functional form of the birth rate. We require analytic solutions of the models or that their potential singularities consist only in removable poles. Such process, known asanalytical modeling, was introduced in [17]. In such work, the required conditions leading to special nonlinearities that discard traditional interaction terms like the law of mass action and logistic terms were established. Examples were exhibited with predation functions where coexistence between the prey population and the predator population was achieved. One possible drawback associated with this new approach is that in order to achieve coexistence, restrictive conditions have to be imposed on the systems. Hopefully, simpler models can be thoroughly analyzed and then built upon to derive increasingly more complex models and this is the basic approach that we will use here.Before attempting to model the vital-rate changes, we need to address the critical lack of age-structured information. To do this, we propose in Section2 a parametrized birth rate which generalizes the one presented in previous works. Then in Section 3, we use a general framework in order to derive predator-prey models with age structure based on the new birth rate. In Section 4, we analyze the models including numerical simulations and interpret the biological implications of the results. Finally, conclusions are given in Section 5.


## 2. Birth Rate Prey Analysis

Population models with age structure are commonly used to analyze the evolution of the dynamics of animal populations with delayed maturation and long life spans, among others. One of the primary building blocks for such population models is the age-specific birth rate. The birth rate is usually the dominant factor in determining the rate of population growth. It depends on the level of fertility, the age structure of the population, and many other factors. The qualitative properties of a birth rate include the fact that it is expected to be small for newborns and old individuals, whereas for young adults it must have the shape of unimodal function with a maximum value for matured adults. In our case we assume that the birth rate depends only on the age of population. Mathematically, a birth rate denoted byβ(ω) is a continuous function defined only for nonnegative ages; that is, ω≥0. Furthermore, it is a positive bounded function with compact support. Since the birth rate for old individuals is 0. That is, β(ω)=0, for ω>A1 and β(ω)<A2 for all w for some constants A1 and A2.There are several functional forms of the functionβ that have been introduced in the literature. For example, in [18], we have introduced an appropriate reproductive rate of the prey population given by (1)βω=β0ωe-αωwith β0>0 and α≥0. This function does not have compact support, but for large values of ω, it is practically zero. Let us recall that if α>0 the birth function (1) behaves in a way that is appropriate for many mammals and if α=0, then β=β0ω has been used for certain fish species. In this case β vanishes for ω equal to zero and it approaches zero for large values of ω. It reaches a global maximum at the age ω=1/α equal to β=(β0/αx)e-1. Here, we propose a parametrized birth rate which includes (1) as a particular case. This new rate is given by(2)βω,m=β0wme-αw,where m is a natural number. Notice that the new rate has mathematically and biologically all the required properties, since it is practically zero for larger values of w. See Figure 1, where we plot the graph of β(ω,m) for some specific values of m. It is important to clarify that the new birth rate, with m, an integer, will allow us to obtain a differential equation system instead of an integrodifferential system as we will discuss in the next section.Figure 1
Graphs of the birth rateβ(ω,m) for different values of m.


## 3. Age-Structured Models

In this section, we develop age-structured models that include our new proposed birth rate function as well as self-limitation of the prey and generic predation interactions. Simplified versions of these new models were previously analyzed in [17, 19]. In those models, a basic assumption was made by considering the birth rate (1). The role of this birth rate is also important because it allows us to obtain the following system of differential equations (instead of an integrodifferential system):(3)dP1dt=-μ0P1-εP12+FP1,P2,P3,P4,dP2dt=-θP2-εP1P2+β0P3,dP3dt=-θP3-εP1P3+FP1,P2,P3,P4,dP4dt=-bP4+GP1,P2,P3,P4,where P1, P2, P3, and P4 represent the total population of the prey, the number of prey newborns that survive cannibalism, the net mortality rate at equilibrium, and the total predator population at time t, respectively. The functions F(P1,P2,P3,P4) and G(P1,P2,P3,P4) are general predation functions and the coefficients are positive constants. There exists a relationship between the functions F and G that for particular cases can be established.Let us remark that this model assumes, as in [20], that newborns are predated by the adult population and only a fraction of them survive this cannibalism per unit of time; that is, (4)FP1,P2,P3,P4=P2tgP1t,P3t,P4t,where g(P1,P3,P4) is the the fraction surviving; hence 0≤g(P1,P3,P4)≤1. The conversion of prey newborns to food is modeled by assuming that the prey death function rate depends on both P1 and P2.Next, we will construct age-structured models based on the birth rate (2). Let us consider that the prey population of age ω at the time t is denoted by ρ(ω,t). The total population is the sum over the ages:(5)P1t=∫0∞‍ρω,tdω.We also consider an age structure by using the McKendrick equation [11]: (6)ρωω,t+ρtω,t+μρω,t=0.Now, let us find differential equations for the total population,P1(t), and the population of newborns, P2(t). In order to do this, we integrate (6) from zero to infinity obtaining: (7)dP1tdt=-μP1P1t+ρ0,t.The newborn population is calculated as P2(t)=∫0∞‍β(w,m)ρ(w,t)dw. Its corresponding differential equation is obtained if we multiply (6) by β0wme-αω and the new birth rate, integrate again, and use integration by parts:(8)dP2tdt=-μP1+αP2t+βom∫0∞‍wm-1e-αwρw,tdw.This integrodifferential equation can be transformed to an ordinary differential equation by defining the variableP3(t) as(9)P3t=∫0∞‍wm-1e-αwρw,tdw.This new variableP3 also satisfies the following differential equations obtained by deriving (9) with respect to t to get(10)dP3tdt=-μ+αP3t+m-1∫0∞‍wm-2e-αwρw,tdw.Again, we get an integrodifferential equation, but the procedure described above can be used by defining(11)Pjt=∫0∞‍wm+2-je-αwρw,tdwforj=2,…,m+2.In general, forj=2,…,m+1 we obtain(12)dPjtdt=-μ+αPjt+m+2-j∫0∞‍wm+1-je-αwρw,tdw,which by definition can be written as (13)dPjtdt=-μ+αPjt+m+2-jPj+1t.Finally forj=m+2 we get(14)dPm+2tdt=-μ+αPm+2t+ρ0,t.Since we have discovered that an essential mechanism to achieve coexistence is the self-limitation of the prey, we will assume thatμ=μ0+ϵP1. Thus, by substituting it in the model we obtain the following prey age-structured model:(15)dP1tdt=-μ0P1-ϵP12+ρ0,t,dP2tdt=-θP2-ϵP1P2+mβ0P3,dP3tdt=-θP3-ϵP1P3+m-1P4,⋮dPm+2tdt=-θPm+2-ϵP1Pm+2+ρ0,t,where μ0, ϵ, η, θ=α+μ0, β0, b, and c are positive parameters.There are two important choices to make in the model. The first one is selecting the incoming of offspring into the prey population,ρ(0,t), and the second is how to postulate the quantitative growth of the predation process. The quantity ρ(0,t) represents the number of newborns that survive cannibalism per unit time. A particular case was given in [8] where it was considered that ρ(0,t)=P2g(P1(t),P2(t)). We make the assumption that ρ(0,t) depends only on Pj(t) for j=1,2,…,m+3 to include different forms of cannibalism. Therefore ρ(0,t)=F(P1,…,Pm+2,Pm+3) for some suitable functions F. In [17, 19], we propose general predation functions suggested by the Painlevé analysis, which we will use in this work. We will denote Pm+3 as the predator population and its interaction given by dPm+3(t)/dt=-bPm+3+G(P1,…,Pm+2,Pm+3). So our final models take the form(16)dP1tdt=-μ0P1-ϵP12+FP1,…,Pm+2,Pm+3,(17)dP2tdt=-θP2-ϵP1P2+mβ0P3,dP3tdt=-θP3-ϵP1P3+m-1P4,⋮dPm+2tdt=-θPm+2-ϵP1Pm+2+FP1,…,Pm+2,Pm+3,dPm+3tdt=-bPm+3+GP1,…,Pm+2,Pm+3.


## 4. Numerical Results

In this section we focus on the numerical analysis of system (17) by choosing appropriate functions F and G, in order to obtain populations coexistence. The function forms of such functions were obtained from [17], they include some known biological models such as the Lotka-Volterra model, the predator-prey model with Holling type predation functions. In our general setting, we do not have, so far, the biological information about the relationship between both functions. We select some examples from [17] in order to illustrate the effect of varying the parameter m. Finding relevant predation functions is a difficult task that depends on the particular system analyzed. It is important to remark that coexistence is established by requiring, not only a stable equilibrium point but also a stable periodic solution. This requirement is of paramount technical difficulty, however commonly found in real life interacting populations.Our goal is to analyze how the birth rate modifies the behavior of system (17), that is, to investigate the effect of varying the power m in the birth rate. In order to do this, we will select some choices of the parameter m and we divide our analysis into several cases.In the first case we chooseF=fP2/(1+Pm+3) considering that cannibalism is only a mean of population control, whereas G=aP2Pm+3 follows the law of mass action. This case has been analyzed in [17] for m=1. In the second case we select F=fP1Pm+3/(1+P2) where cannibalism is primary source of food and G=aPm+3/(P2+σ1), in the third case F=fP2/(1+Pm+3) and G=aPm+3/(1+P3), and finally, we take in the fourth case F=fP1P2/(σ2+Pm+3) and G=aP22/(σ3+Pm+3). All of these combinations have presented periodic branches of solutions for the case m=1.


### 4.1. Case 1 (F=fP2/(1+Pm+3) and G=aP2Pm+3)

This combination of functionsG and F corresponds to a newborn predation interaction of the Lotka-Volterra type and a saturation term, respectively. The function F was first obtained in [18] with the following approximation assuming that the predator population, Pm+3, is small:(18)FP2,Pm+3=P2t-μ0P2tPm+3t=P2t1+μ0Pm+3t+OPm+32t.We take F as the following proportional truncated form assuming that μ0=1 and coefficient f is added to weigh the contribution of F in the model: (19)FP2,Pm+3=fP2t1+Pm+3t.The choice of functionG, the Lotka-Volterra predation term, was made by considering the simplest case for the predation interaction. System (17) with this selection of functions F and G with m=1 and fixed coefficients was proposed as a viable solution for coexistence between the species. Thus, it is an important case to analyze by modifying the prey birth rate.We will fix every parameter exceptα, with values μ=0.77125, ε=0.01, f=3.5416, β0=0.46875  b=0.03, and a=0.10. As we vary α, we find that there are eleven solution branches for every value of α and every value of m. Of those branches, only five are of biological interest. Four of them correspond to equilibrium branches and the other one is a periodic branch. The stability of equilibria varies with the parameter α. However, what is of paramount importance is the fact that for m=1 such periodic branch is unstable but for m>1 it is stable. Following the stable periodic branch solution (m>1), the period increases up to infinity. In biological terms these results are important since by varying the birth rate we are able to obtain coexistence from an unstable system. These facts can be seen more clearly from Figure 2 where we plot the parameter α against the L2 norm of the solution. For stationary solutions we use the Euclidean norm and for periodic solutions we use the norm(20)1T∫0T‍∑j=1m+3‍Pj2tdt,where T is the period of the solution. Branches of stable periodic solutions are indicated by black curves and the stationary branches by colored curves. A dotted curve means that the solution is unstable for those values of the parameter.Figure 2
Case 1. (a)m=2, (b) m=5.
(a)
(b)


### 4.2. Case 2 (F=fP1Pm+3/(1+P2) and G=aPm+3/(P2+σ1))

In this case we study the case of a type of social prey which becomes aggressive towards the predators when their number is large enough but becomes susceptible for small numbers; for example, we can consider the African buffalo (Syncerus caffer) as the prey. Thus we choose G proportional to Pm+3/(P2+σ1). For F we select the truncated form F=fP1Pm+3/(1+P2) analogous to (19). This time the calculation of the continuation of equilibrium points and periodic solutions become more complicated. We again fix the values of all the parameters and vary α. The values are the same as in Case 1, except that ε=0 meaning that self-limitation is not considered. We find that there is only one branch of equilibrium solutions and one branch of periodic solutions for every value of α and only for m=1,2 and for every particular initial condition. The periodic branch is stable and its period also increases to infinity. In the case of m=1 the periodic branch changes stability only for a small parameter regime. Notice that even for the case m=1, coexistence is established. In Figure 3 we can observe these facts in detail. What is remarkable in this case is that self-limitation is not necessary to achieve coexistence.Figure 3
Case 2. (a)m=1, (b) m=2.
(a)
(b)


### 4.3. Case 3 (F=fP2/(1+Pm+3) and G=aP2Pm+3/(1+P3))

In this case we chooseG proportional to P2Pm+3/(1+P3). For F we select the truncated form F=fP2/(1+Pm+3) analogous to (19). The values of the constants are the same as in Case 1. As in Case 2, there is only one branch of equilibrium solutions and one branch of periodic solutions at least for the cases m=1,2,3. Bigger values of m were not analyzed. The periodic branch is stable and its period increases faster to infinity than those for Cases 1 and 2. In this case m acts as a perturbation parameter for the Hopf bifurcation points (see Figure 4).Figure 4
Case 3. (a)m=1, (b) m=3.
(a)
(b)


### 4.4. Case 4 (F=fP1P2/(σ2+Pm+3) and G=aP22/(σ3+Pm+3))

This case is presented only for theoretical purposes in order to illustrate the possibilities in the selection of the predator interaction and the fraction of newborns that survive predation. The qualitative behavior is the same as in Case 3 and the only difference is the variation of the equilibrium branch; in this case we choosef as the variation parameter (see Figure 5).Figure 5
Case 4. (a)m=1, (b) m=3.
(a)
(b)


## 4.1. Case 1 (F=fP2/(1+Pm+3) and G=aP2Pm+3)

This combination of functionsG and F corresponds to a newborn predation interaction of the Lotka-Volterra type and a saturation term, respectively. The function F was first obtained in [18] with the following approximation assuming that the predator population, Pm+3, is small:(18)FP2,Pm+3=P2t-μ0P2tPm+3t=P2t1+μ0Pm+3t+OPm+32t.We take F as the following proportional truncated form assuming that μ0=1 and coefficient f is added to weigh the contribution of F in the model: (19)FP2,Pm+3=fP2t1+Pm+3t.The choice of functionG, the Lotka-Volterra predation term, was made by considering the simplest case for the predation interaction. System (17) with this selection of functions F and G with m=1 and fixed coefficients was proposed as a viable solution for coexistence between the species. Thus, it is an important case to analyze by modifying the prey birth rate.We will fix every parameter exceptα, with values μ=0.77125, ε=0.01, f=3.5416, β0=0.46875  b=0.03, and a=0.10. As we vary α, we find that there are eleven solution branches for every value of α and every value of m. Of those branches, only five are of biological interest. Four of them correspond to equilibrium branches and the other one is a periodic branch. The stability of equilibria varies with the parameter α. However, what is of paramount importance is the fact that for m=1 such periodic branch is unstable but for m>1 it is stable. Following the stable periodic branch solution (m>1), the period increases up to infinity. In biological terms these results are important since by varying the birth rate we are able to obtain coexistence from an unstable system. These facts can be seen more clearly from Figure 2 where we plot the parameter α against the L2 norm of the solution. For stationary solutions we use the Euclidean norm and for periodic solutions we use the norm(20)1T∫0T‍∑j=1m+3‍Pj2tdt,where T is the period of the solution. Branches of stable periodic solutions are indicated by black curves and the stationary branches by colored curves. A dotted curve means that the solution is unstable for those values of the parameter.Figure 2
Case 1. (a)m=2, (b) m=5.
(a)
(b)


## 4.2. Case 2 (F=fP1Pm+3/(1+P2) and G=aPm+3/(P2+σ1))

In this case we study the case of a type of social prey which becomes aggressive towards the predators when their number is large enough but becomes susceptible for small numbers; for example, we can consider the African buffalo (Syncerus caffer) as the prey. Thus we choose G proportional to Pm+3/(P2+σ1). For F we select the truncated form F=fP1Pm+3/(1+P2) analogous to (19). This time the calculation of the continuation of equilibrium points and periodic solutions become more complicated. We again fix the values of all the parameters and vary α. The values are the same as in Case 1, except that ε=0 meaning that self-limitation is not considered. We find that there is only one branch of equilibrium solutions and one branch of periodic solutions for every value of α and only for m=1,2 and for every particular initial condition. The periodic branch is stable and its period also increases to infinity. In the case of m=1 the periodic branch changes stability only for a small parameter regime. Notice that even for the case m=1, coexistence is established. In Figure 3 we can observe these facts in detail. What is remarkable in this case is that self-limitation is not necessary to achieve coexistence.Figure 3
Case 2. (a)m=1, (b) m=2.
(a)
(b)


## 4.3. Case 3 (F=fP2/(1+Pm+3) and G=aP2Pm+3/(1+P3))

In this case we chooseG proportional to P2Pm+3/(1+P3). For F we select the truncated form F=fP2/(1+Pm+3) analogous to (19). The values of the constants are the same as in Case 1. As in Case 2, there is only one branch of equilibrium solutions and one branch of periodic solutions at least for the cases m=1,2,3. Bigger values of m were not analyzed. The periodic branch is stable and its period increases faster to infinity than those for Cases 1 and 2. In this case m acts as a perturbation parameter for the Hopf bifurcation points (see Figure 4).Figure 4
Case 3. (a)m=1, (b) m=3.
(a)
(b)


## 4.4. Case 4 (F=fP1P2/(σ2+Pm+3) and G=aP22/(σ3+Pm+3))

This case is presented only for theoretical purposes in order to illustrate the possibilities in the selection of the predator interaction and the fraction of newborns that survive predation. The qualitative behavior is the same as in Case 3 and the only difference is the variation of the equilibrium branch; in this case we choosef as the variation parameter (see Figure 5).Figure 5
Case 4. (a)m=1, (b) m=3.
(a)
(b)


## 5. Conclusions

The main contributions of this work are that (a) we have presented a new series of age-structured models in order to show the importance of the birth rate for coexistence, (b) we have shown that by changing the rate of birth there exist periodic solutions for such models, and (c) in order to establish coexistence, we sought stable equlibrium branches and more importantly stable periodic solution branches. We also find out that there are complex models where the structure of their bifurcation diagrams were not significantly changed from some simpler models. All the numerical studies were carried out using different predation functions obtained by analytical techniques from the literature. We conclude that the use of more general rates of birth is a necessary mechanism to achieve realistic models. As a practical applicability of this work, it is worth mentioning that there exist a variety of biological control methods to prevent the spread of pests which can be analyzed with our models. The idea is to control the population by finding stable periodic solutions of the populations involved.



---

*Source: 241312-2015-08-05.xml*