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

and Applied Analysis 3 is modeled by assuming that the prey death function rate depends on both P 1 and P 2 . 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:


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][2][3][4][5][6][7][8][9][10]. One relevant feature in these models is the predator-prey coexistence, which is of vital importance to obtain realistic models of predatorprey 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 Mac-Camy [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 as analytical 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.

2
Abstract and Applied Analysis 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 Section 2 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.

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 > 1 and ( ) < 2 for all for some constants 1 and 2 .
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 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 Here, we propose a parametrized birth rate which includes (1) as a particular case. This new rate is given by where 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 . See Figure 1, where we plot the graph of ( , ) for some specific values of . It is important to clarify that the new birth rate, with , an integer, will allow us to obtain a differential equation system instead of an integrodifferential system as we will discuss in the next section.

Age-Structured Models
In this section, we develop age-structured models that include our new proposed birth rate function as well as selflimitation 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): where 1 , 2 , 3 , and 4 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 , respectively. The functions ( 1 , 2 , 3 , 4 ) and ( 1 , 2 , 3 , 4 ) are general predation functions and the coefficients are positive constants. There exists a relationship between the functions and 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, The conversion of prey newborns to food is modeled by assuming that the prey death function rate depends on both 1 and 2 . 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 is denoted by ( , ). The total population is the sum over the ages: We also consider an age structure by using the McKendrick equation [11]: Now, let us find differential equations for the total population, 1 ( ), and the population of newborns, 2 ( ). In order to do this, we integrate (6) from zero to infinity obtaining: The newborn population is calculated as 2 ( ) = ∫ ∞ 0 ( , ) ( , ) . Its corresponding differential equation is obtained if we multiply (6) by 0 − and the new birth rate, integrate again, and use integration by parts: This integrodifferential equation can be transformed to an ordinary differential equation by defining the variable 3 ( ) as This new variable 3 also satisfies the following differential equations obtained by deriving (9) with respect to to get Again, we get an integrodifferential equation, but the procedure described above can be used by defining In general, for = 2, . . . , + 1 we obtain which by definition can be written as Finally for = + 2 we get Since we have discovered that an essential mechanism to achieve coexistence is the self-limitation of the prey, we will assume that = 0 + 1 . Thus, by substituting it in the model we obtain the following prey age-structured model: where 0 , , , = + 0 , 0 , , and 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, ), and the second is how to postulate the quantitative growth of the predation process. The quantity (0, ) represents the number of newborns that survive cannibalism per unit time. A particular case was given in [8] where it was considered that (0, ) = 2 ( 1 ( ), 2 ( )). We make the assumption that (0, ) depends only on ( ) for = 1, 2, . . . , + 3 to include different forms of cannibalism. Therefore (0, ) = ( 1 , . . . , +2 , +3 ) for some suitable functions . In [17,19], we propose general predation functions suggested by the Painlevé analysis, which we will use in this work. We will denote +3 as the predator population and its interaction given by +3 ( )/ = − +3 + ( 1 , . . . , +2 , +3 ). So our final models take the form

Numerical Results
In this section we focus on the numerical analysis of system (17) by choosing appropriate functions and , 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 . 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 in the birth rate. In order to do this, we will select some choices of the parameter and we divide our analysis into several cases.

Case 1 (
). This combination of functions and corresponds to a newborn predation interaction of the Lotka-Volterra type and a saturation term, respectively. The function was first obtained in [18] with the following approximation assuming that the predator population, +3 , is small: We take as the following proportional truncated form assuming that 0 = 1 and coefficient is added to weigh the contribution of in the model: The choice of function , the Lotka-Volterra predation term, was made by considering the simplest case for the predation interaction. System (17) with this selection of functions and with = 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, = 3.5416, 0 = 0.46875 = 0.03, and = 0.10. As we vary , we find that there are eleven solution branches for every value of and every value of . 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 = 1 such periodic branch is unstable but for > 1 it is stable. Following the stable periodic branch solution ( > 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 2 Abstract and Applied Analysis norm of the solution. For stationary solutions we use the Euclidean norm and for periodic solutions we use the norm where 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.

4.2.
Case 2 ( = 1 +3 /(1 + 2 ) and = +3 /( 2 + 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 proportional to +3 /( 2 + 1 ). For we select the truncated form = 1 +3 /(1 + 2 ) 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 = 1, 2 and for every particular initial condition. The periodic branch is stable and its period also increases to infinity. In the case of = 1 the periodic branch changes stability only for a small parameter regime. Notice that even for the case = 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.
For we select the truncated form = 2 /(1 + +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 = 1, 2, 3. Bigger values of 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 acts as a perturbation parameter for the Hopf bifurcation points (see Figure 4).

4.4.
Case 4 ( = 1 2 /( 2 + +3 ) and = 2 2 /( 3 + +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 choose as the variation parameter (see Figure 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.