
MBW:Three Species Food Chain ModelingFrom MathBioWilliam Dowling, Toni Klopfenstein, Lucas Koepke ContentsSummarySongjuan Lv and Min Zhao examine a threespecies food chain with a Hollingtype functional response in the paper The Dynamic Complexity of a Three Species Food Chain Model ^{[1]}. By analytically examining and numerically simulating the system, the authors analyze the food chain’s chaotic, stable and periodic dynamics. The authors examine the four nonnegative equilibrium points for the system, and perform a local stability analysis for each equilibrium point. The longterm behavior of the system is also examined using numerical analysis. It is also shown how the system can be pushed into chaos or into a steady state by varying certain parameters, such as the prey carrying capacity or the prey growth rate. An examination of the largest Lyapunov exponent is also included. BackgroundThe three species food chain model is an extension of the general two species model given by Lotka and Volterra. Until recently, it was believed that only steady states and limiting cycles were exhibited in those types of systems. However, researchers have recently discovered the appearance of chaos and cycles in predatorprey species models with nonlinear responses, including BeddingtonDeAngelis ^{[2]} and Hollingtype responses. This paper deals primarily with the Holling functional response. C.S. Holling^{[3]} ^{[4]}discovered this response to describe the consumption rate of prey and the relationship with the prey population density. There are three different characteristic responses included in the Holling functional response family.^{[5]} In Type I responses, prey consumption grows linearly until reaching an equilibrium point at which consumption maintains a constant value. The consumption rate increases at the rate of predator attack rates. These are usually found in a herbiovoreoplant or invertebrate predatorprey food chains^{[6]} In Type II responses, prey consumption grows with an increase in prey density until a constant consumption level is reached, regardless of any other increases in prey density. This is the most commonly observed functional response in nature, and is modeled by Holling's Disk Equation, ^{[7]} shown below, (This asymptotic response is more commonly seen in terms of N, as in APPM4390:HostParasitoid_Models#Type_II:_Asymptotic_in_N, and is also described in discrete time by the Holling Model  APPM4390:HostParasitoid_Models#Holling_Model) Type III responses have an initial increase in prey consumption which then levels out to a constant consumption level. The system examined in this paper involves a threespecies trophic system, in which there are two competing predators (one intermediate and one top predator) and one prey species, and has a Holling TypeII functional response. Examples of this type of system include mousesnakeowl, wormrobinfalcon, and vegetationharelynx. ^{[8]} Data taken from wild food chains is often incomplete and so it can be difficult to verify models. However in certain cases it is possible to reconstruct the states of all the species in a three species model. For more about food chain modeling see APPM4390:Robustness Analysis of an ObserverBased Controller in a Food Web Mathematical ModelThe general system of equations used to represent a three species food chain is: given by Varriale and Gomes.^{[9]} The parameters are defined as: In this equation, the functions F_{1} and F_{2} are representations of the interactions between species, which represent a Hollingtype II functional response, and are given by the equations below. If the original system is modified to include this Hollingtype functional response, the following model is achieved. The new parameters are: The model can then be transformed using: and the dimensionless model is created: with the new parameters: Equilibria AnalysisStability analysis of equilibria is performed on the nondimensionalized form of the system given above. There are four possible nonnegative equilibrium points. They correspond to the following biological interpretations: (1) all species extinct; (2) prey population at carrying capacity in absence of both predators; (3) top predator absent and intermediate predator survives on its prey; (4) all species present. The equilibrium points (1) and (2) corresponding to (0,0,0) and (k0,0,0) always exist. The Jacobian matrices of both are given below: The eigenvalues of both are given by the diagonal entries in the corresponding Jacobian matrix. Both have one nonnegative and two negative eigenvalues, from which we determine that both equilibria are saddle points. (Gakkhar et al similarly found saddle point equilibria when looking at delayed preypredator models: APPM4390:AliHussain#Result) The equilibrium point (x*,y*,0) corresponding to (3) where x* and y* are defined as exists only when the following conditions are met: In order to study the stability of (3), we can study the following two dimensional subsystem for which (3) is an equilibrium point: The Jacobian of this system at (x*,y*) is given by: This equilibrium is asymptotically stable if both coefficients of the corresponding characteristic equation are greater than zero. The coefficients for the characteristic equation are given by The conditions given for the existence of an equilibrium at (x*,y*,0) force the second coefficient to be greater than zero. Below are two graphs showing the change in the first coefficient with respect to r0 and k0: Note that for the range of values tested, the value of the first coefficient never surpasses zero except for k0 < 1.3 (roughly). As a result, we can determine that the equilibrium point is not asymptotically stable for all values tested except for k0 < 1.3. However, the condition 0 < x* < k0 < k1 must still hold for the equilibrium point to exist (for the parameter used here, x* = 1/3). The fourth equilibrium exists for positive x, y, and z values if and only if there exists a positive solution to the following system of equations (derived directly from the nondimensionalized system of differential equations being studied): The solution is given by where The Jacobian at (4) is given by where The coefficients of the corresponding characteristic equation are given by The equilibrium point (4) is locally stable if the following condition is met: Below, the quantities and are graphed for a range of r0 and k0 values (the median possible x value is used): Note that all must be greater than zero in order for the equilibrium point to be asymptotically stable. This does not occur for any of the r0 or k0 values tested, thus we can determine that the equilibrium (4) is not asymptotically stable for any of these values. Numerical AnalysisTo study the behavior of the nondimensionalized system, we will examine how varying one of the input parameters, while keeping the others fixed, affects the system. As an example, we look at the parameter r_{0}, which is the growth rate of species x, and how it determines the characteristics of the top predator z. To do this, we will use a Poincare map of the local maxima for the top predator z. The fixed parameter values are: with initial condition x_{0} = 3, y_{0} = 0.1 and z_{0} = 0.2 Poincare mapThe Poincare map is a useful tool in the analysis of dynamical systems, and is defined as follows. Given an ndimensional system, let S be an n1 dimensional surface such that all flow through S is entirely perpendicular to it. That is, there are no trajectories, starting on S, with flow parallel to S. Then the Poincare map, call it P, is defined as a mapping from S to S, such that x_{k+1} = P (x_{k}) Then a closed orbit of the original system will be a fixed point of P, such that after some time t, P (x*) = x*. In practice, it is nearly impossible to explicitly determine this mapping P, however using numerical techniques the Poincare map provides a useful framework to determine when cycles exist in a dynamical system. Specifically, the continuous solution to the original equation will have to be discretized in some way to form a sequence a values that can be analyzed for points that fit the definition for a fixed point for P. A natural way to do this is to look at the sequence of values corresponding to the local maxima, and any values that repeat will be fixed points of P. Numerical resultsTo use this technique on the nondimensionalized system given above, we will let r_{0}, which is the growth rate of the prey x, vary between 1.5 and 9, and then analyze the changes in behavior produced in the population of the top predator z. A sequence of the local maxima of the time series for z is generated, and then any repeated values in this sequence show a cycle exists in the original system. To begin with, for r_{0} = 1.5, the plot of z with respect to time looks like and highlighting the maxima with red dots looks like Starting from the left side of the picture, we let x_{1} be the first maxima, and then x_{2} = P(x_{1}) be the second, and so forth. The sequence of values generated for the local maxima is then {1.416, 1.076, 1.207, 1.199, 1.201, 1.201, 1.2, 1.201, 1.201, 1.201, 1.201, 1.201} And since 1.201 is repeated, and is in fact the only repeated value, we now have a single fixed point for P, which corresponds to a period1 cycle in the original system.
and highlighting the maxima with red dots looks like Then the sequence of values generated for the local maxima is {1.909, 0.957, 2.091, 0.778, 2.082, 0.785, 2.081, 0.786, 2.081, 0.787, 2.081, 0.787, 2.081, 0.787, 2.081} which has two fixed points for P, which correspond to a period2 cycle in the original system.
As r_{0} is further increased, we see period doubling, first to a period4 cycle, and then chaos until r_{0} = 6.5 where a period3 cycle emerges, only to undergo period doubling again, leading to chaos. ConclusionsAs shown by the stability analysis, we determined that separately varying two of the parameters, r_{0} and k_{0}, within the set limits would never allow an asymptotically stable solution. This was first shown using the RouthHurwitz criteria for the coefficients in the characteristic equation, and since one of them was always less than zero, no asymptotically stable solution could exist. This was verified with the numerical results for r_{0} between 1.5 and 9, since there was always at least one fixed point for the Poincare map, and hence there was always at least one cycle in the original system. This periodic nature of the solutions to our model for three interacting species fits our intuition, since we would not expect to have an exactly constant number of any one species in the real world. Since we were able to reproduce and verify the results found in the paper we reviewed, using Matlab and Mathematica, we conclude that a justifiable and accurate model was used. The methods for both stability analysis and numerical modeling are techniques commonly applied to dynamical systems, and the results from both of these methods showed the same conclusion. The occurrence of bifurcations and chaos are also possible with a three species model, and both of these cases were observed using numerical modeling. As a result, studying a three species food chain yields complex dynamics and interesting behavior, lending itself well to mathematical insight, while still staying grounded in reality. Project categorizationMathematicsThis model uses a system of ordinary differential equations to model the population dynamics of three interacting species. Stability analysis of the found equilibrium points is performed using Jacobian matrices and eigenvalues. Model typeThe system of ODEs are modeling the change in each of the different species populations with respect to time. Biological systemThe biological system studied is the interaction of three competing species. Related TopicsAnother example of useful population modeling techniques would be the use of Trophic pyramids which don't count individuals, but rather net biomass or production rates. Recent CitationThe above article was recently cited in Min Zhao and Limin Zhang’s 2009 paper titled “Permanence and chaos in a host–parasitoid model with prolonged diapause for the host”. The paper discusses the effects of a prolonged dormant period by the parasitoid, and how it alters the behavior of the host parasitoid relationship. The relationship model is described by a pair of coupled difference equations. The permanence of the system is investigated through bifurcation analysis and using the largest Lyapunov exponents while altering the diapause rate d. It was found that when 0.281 < d < 0.4879, the Lyapunov exponents are negative, which means that there exist a stable coexistence state between the host and parasitoid. The article was cited when looking at the complexities of the observed system. The parasitoid host model has similar complex dynamic patterns found in the bifurcation analysis as the three – species foodchain model. Sources
