Editorial: Dynamics of Deterministic Models of Biological Systems

1.   Introduction

Trophic interactions are the foundation of biology and allow us to deepen our understanding of the interplay between life forms. Predation is probably the most well-known of these interactions, describing how two populations interact as hunters (predators) and food source (prey). When predation occurs among organisms of the same species, it is called cannibalism. Cannibalism is of particular interest because of its paradoxical nature and prevalence among 1300 species ^[1] including several species of amphibian larvae (tadpoles) ^[2], fish ^[1,3,4], and birds of prey ^[5]. Cannibalism can occur at multiple trophic levels within a single biome and can be isolated to a single life history stage ^[1]. For instance, tadpoles and hatchling birds at their juvenile stages of life will cannibalize conspecifics under certain conditions, but as they mature, this behavior is rarely seen ^[2,5]. Many biological studies have investigated the motivations behind cannibalism and found various gains to the organism ^[1], some of which includes life history benefits by working to increase an individual's reproductive fitness. For example, when a succeeding male lion takes over a female pride, it may cannibalize the unrelated young in order to copulate and sire offspring sooner, increasing its reproduction success ^[6]. Cannibalism can also increase an individuals chances of survival, especially in nutrient-depleted environments. In some species of tadpoles, consuming conspecifics increases the tadpole's rate of growth, allowing the tadpole to mature at a faster rate and leave the ephemeral pool before it dries out ^[2]. Studies in ^[7] also showed that cannibalism contributed to the reduction in growth rate of Parrotfish and Pterois Volitans population. Despite this, opportunistic cannibalism offers an additional food source with little risk of injury to the cannibal ^[1]. For example, Atlantic cod, (Gadus marhua), are found to cannibalize within the larval stage since relative size differences between batches of hatchlings allows larval cod to successfully attack, handle, and consume smaller fellow young ^[4]. In addition, cannibalism can also decrease competition. In many birds of prey, such as the common buzzard, (Buteo buteo), older hatchlings will cannibalize younger chicks in order to reduce competition and increase parental care ^[5]. Dragonfly larvae exhibited density-dependent cannibalism, consuming smaller conspecifics more frequently in high density environments, aiding in the regulation of the population and emergence synchrony ^[8]. In recent years, novel incidents of cannibalism have been noted in glaucous-winged gulls (Larus glaucescens) ^[9] and polar bears (Ursus maritimus) ^[10] due to global warming destroying habitats and increasing competition between conspecifics. One well known effect of cannibalism is the so called "life-boat mechanism" where cannibalistic predators are capable of surviving during food shortage periods. In such circumstances, non-cannibalistic predators are not able to survive ^[11,12].

Mathematical models have been used to gain insights into the ramifications of this macabre behavior. In ^[13], cannibalism was modeled in both predator and prey populations in a two species Holling-Tanner predator-prey system. The model is given by

where $u$ and $v$ represent the prey and predator population respectively. The term $c \left(\dfrac{u^2}{u+d}\right)$ models the functional response of the cannibalistic prey and $c_1$ the energy gained from cannibalism. Also, $\rho v$ is the cannibalism term in the predator population with a rate of cannibalism $\rho$. We refer the reader to ^[13] for detailed description of the parameters in the model. Findings showed oscillatory dynamics in the system stabilized when cannibalism was present in both the predator and prey populations, but not when acting separately. Also, results from ^[14] showed that the presence of cannibalism in the prey population had a stabilizing effect when increased above a threshold. Magnusson studied cannibalism in a structured three-species predator-prey system and found that cannibalism had a destabilization effect ^[15]. Similar results were obtained by Kaewmanee and Tang ^[16]. Cushing also showed that cannibalism has both a stabilizing and destabilizing effect depending on the choice of parameters ^[11]. In ^[17], a Lotka-Volterra predator-prey model was studied and results showed that cannibalism had both a positive and negative impact on the stability of the system. Kohlmeier and Ebenhöh in ^[14] also showed that, cannibalism in predators is a stabilizing mechanism in a predator-prey age structured model. Spatially explicit models have also been studied to explore the impacts of cannibalism ^[18,19,20]. Cannibalism was found to have a destabilizing effect in a spatially explicit three-species age structured predator-prey system, causing the emergence of Turing patterns ^[19]. Spatial patterning results were obtained in ^[20] when cannibalism was modeled in the predator population. There has been reported cannibalism among tadpoles in a specie of treefrog (Osteopilus septentrionalis) ^[21]. In the concluding work of Aladeen in ^[13], there are no mathematical models that explore cannibalism completely competition motivated when there is a high density of tadpoles constrained by small space or less resources in a pond. As such, a stage-structured model where cannibalism is modeled in the juvenile prey best encapsulates this situation. Thus, the focus of this paper is to propose a theoretical stage-structured system of ordinary differential equations to model the relationship between predators, adult prey, and cannibalistic juvenile prey.

In the current manuscript, we report the following:

● Rich dynamical structure is revealed via Bogdanov-Takens and cusp bifurcations respectively using numerical examples. See Figures 1 and 2.

● Cannibalism has a stabilization effect. Increasing the cannibalism rate of attack changes the interior equilibrium from an unstable to a stable steady state as seen in Figures 3 and 4.

● Cannibalism has a destabilization effect. Increasing/decreasing the cannibalism rate of attack changes the interior equilibrium from a stable to an unstable steady state as seen in Figure 5.

● In the presence of cannibalism, the system experiences a Hopf bifurcation via Theorem 3.1 and a saddle-node bifurcation via Theorem 3.2.

● Limit cycle dynamics do not exist in prey populations in the absence of predators via Theorem 5.1.

● The presence of cannibalism can induce coexistence of both prey and predator populations. The absence of it leads to extinction of both populations. See Figure 7.

The rest of the paper is organized as follows: Section 2 deals with the formulation and mathematical analysis of our model. We discuss local bifurcation analysis in Section 3. In Section 4, we perform numerical experiments to corroborate our analytical results. We study the dynamics of a subsystem of the formulated model when there are no predators in Section 5. Finally in Section 6, we present a discussion and conclusion of our results.

2.   Model formulation and analysis

We consider a stage-structured predator-prey model where the prey population is categorized into two classes, namely, the juvenile prey and adult prey. We only model cannibalism in the juvenile prey population. The juvenile prey, adult prey and predator population are represented by the state variables $x$, $y$ and $z$ respectively at any instant time $t$. Our model is based on the following ecological assumptions:

● Juvenile prey die as a result of cannibalism and natural death.

● Juvenile prey cannot reproduce and their rate of growth is dependent on the adult prey. Also the growth of the juvenile prey is inhibited by the adult prey through intraspecific competition.

● Predators attack only adult prey.

● When predators are in abundance, density dependent effects such as diseases, overcrowding and intraspecific competition causes quadratic mortality rates ^[22,23,24]. We include a quadratic mortality rate for the predator population in our model.

● Holling Type Ⅱ functional response is used in describing the interaction between predators and adult prey. It is also used in modeling cannibalism in juvenile prey at a rate $c$.

● The energy gained from cannibalism in the juvenile prey is modeled via a $c_1 x$ term where $c_1 < c$.

● The conversion rate of food biomass to predator is always less than the rate at which predators attack and kill adult prey.

2.1.   Model equation

The nonlinear system of ordinary differential equations satisfying our assumptions is given by

with positive initial conditions $x(0) = x_0$, $y(0) = y_0$ and $z(0) = z_0$. All parameters used are assumed positive and their descriptions are given in Table 1.

2.2.   Mathematical analysis of model

The positivity and boundedness of solutions to system (2.1) are important properties for biological meaningfulness. The positivity of solutions implies that the population will continue to thrive and, since resources are limited the population cannot grow beyond certain limits and hence are bounded. We discuss these properties in this subsection.

2.2.1.   Positivity of solutions

We recap the following result which guarantees the positivity of solutions from ^[25,26].

Lemma 2.1. Consider the following system of ODEs:

Non-negativity of solutions is preserved with time, that is,

if and only if

we have

2.2.2.   Boundedness

As stated earlier, for system (2.1) to be biologically meaningful, all solutions initiating from $\mathbb{R}_+^3$ must be bounded. We state the following lemma:

Lemma 2.2. Solutions to system (2.1) are bounded.

Proof. Let us consider the function $W(t) = x(t)+y(t)+z(t)$. Then,

where $\eta = {\rm{min}}((\delta -r), (\mu-c_1))$ and $r\leq \delta, c_1 \leq \mu.$ Then adding $\eta z$ to both sides yields

This implies that

by standard theory in differential inequality. As $t \rightarrow \infty$,

Therefore by Lemma 2.1 and (2.2), all solutions of (2.1) with initial conditions $x(0) > 0, \; y(0) > 0$ and $z(0) > 0$ are bounded in the region

Hence proof.

Remark 1. We remark that the conditions in the proof of Lemma 2.2 are sufficient for system (2.1) to be bounded.

2.3.   Model without cannibalism $(c=c_1=0)$

We first study system (2.1) when there in no cannibalism in the juvenile prey population. This will enable us to understand the impacts of cannibalism when it is present. Therefore, system (2.1) reduces to

2.3.1.   Equilibria and local stability analysis

The system (2.3) has the following biologically feasible non-negative equilibrium points obtained by solving the system $\dfrac{dx}{dt} = 0, \dfrac{dy}{dt} = 0\; \mbox{and}\; \dfrac{dz}{dt} = 0$.

(a) $E_0 = (0, 0, 0)$, which corresponds to the extinction of all populations,

(b) $E_1 = (x_1^*, y_1^*, 0)$, where

(c) $E_2 = (x_2^*, y_2^*, z_2^*),$

where

and $y_2^*$ is a positive real root of the cubic equation

where

The predator free equilibrium $E_1$ exists provided $rm > \delta(\mu +m)$. We provide the conditions for the existence of the interior equilibrium point $E_2$ using the Descartes rule of signs for Eq (2.4).

$({\rm{i}})$ If $m(r-\delta)-\mu \leq 0$, then we cannot find a positive root for Eq (2.4) and therefore $E_2$ does not exist.

$({\rm{ii}})$ If $m(r-\delta)-\mu > 0$, $\kappa_2 < 0$ and $\kappa_3 < 0$ or

$({\rm{iii}})$ $m(r-\delta)-\mu > 0$, $\kappa_2 > 0$ and $\kappa_3 > 0$ or

$({\rm{iv}})$ $m(r-\delta)-\mu > 0$, $\kappa_2 < 0$ and $\kappa_3 > 0$, then there is a unique positive root.

$({\rm{ii}})$ If $m(r-\delta)-\mu > 0$, $\kappa_2 > 0$ and $\kappa_3 < 0$, then there are three positive roots.

Next, we calculate the Jacobian of system (2.3) and obtain

We state the following results pertaining to the local and global stability of the equilibria of system (2.3):

Theorem 2.3. The following statements hold:

(i) The trivial equilibrium point $E_0$ is globally stable if $r < \delta$ ^[27].

(ii) The predator free equilibrium $E_1$ is locally unstable if $rm > \delta(\mu +m)$.

(iii) The positive interior equilibrium $E_2 = \left(x_2^*, y^*_2, z_2^*\right)$ is locally stable provided that the Routh-Hurwitz stability criteria are satisfied.

Proof. We begin proof of global stability for $E_0$ by considering the Lyapunov function $\Xi(x, y, z) = x+y+z$ which satisfies $\Xi(E_0) = 0$ and $\Xi(x, y, z) > 0$ if $(x, y, z)\neq 0$. We also suppose that $r < \delta$. Computing the derivative of $\Xi$ with respect to $t$ yields

since $r < \delta$. Hence the equilibrium point $E_0$ is globally stable.

Next, we evaluate the Jacobian in Eq (2.5) at $E_1$ and obtain

It's associated characteristic equation is

Clearly, since the real part of at least one of the eigenvalues $\lambda_1 = \dfrac{\alpha _1 y^*}{d_1+y^*}$ is positive, we conclude that the predator free equilibrium $E_1$ is locally unstable.

We also note that the characteristic equation of the Jacobian $J^*$ evaluated at $E_2$ is given by

where

The positive interior equilibrium $E_2 = \left(x_2^*, y^*_2, z_2^*\right)$ is locally stable provided that $\beta_1 > 0$, $\beta_2 > 0$, $\beta_3 > 0$ and $\beta_1 \beta_2-\beta_3 > 0$ which satisfies the Routh-Hurwitz stability criteria.

2.4.   Model with Cannibalism $(c > c_1 > 0)$

2.4.1.   Equilibria and local stability analysis

We study the dynamics of system (2.1) when cannibalism is introduced at some rate $c$. The system possesses the following non-negative equilibria. These are

$({\rm{a}})$ $E^*_0 = (0, 0, 0)$,

$({\rm{b}})$ $E^*_1 = (x_1^*, y_1^*, 0),$

where

and $x_1^*$ is a positive real root of the quadratic equation

where

Let $\Gamma = m\left(r-\delta\right)-\delta (\mu -c_1)$. We provide the conditions for the existence of a positive root for Eq (2.9) using the Descartes rule of signs.

$({\rm{i}})$ If $\Gamma \leq 0$, then Eq (2.9) has no positive root.

$({\rm{ii}})$ If $\Gamma > 0$ and either $\dfrac{rmd}{\delta K}+c < \Gamma$ or $\dfrac{rmd}{\delta K}+c > \Gamma$, then there exists a positive root for Eq (2.9).

(c) $E^*_2 = (x_2^*, y_2^*, z_2^*),$

where

and $y_2^*$ is a positive real root of the sextic equation

The constants $\eta_5$, $\eta_4$, $\eta_3$, $\eta_2$, $\eta_1$ and $\eta_0$ are found in the Appendix. We are guaranteed of a positive root for Eq (2.10) if any $\eta_k > 0$ for $k = 0, ..., 5$ by the Descartes rule of signs.

Calculating the Jacobian of system (2.1), we get

We state the following results with regards to the local stability of the equilibria of system (2.1) when cannibalism is present:

Theorem 2.4. The following statements hold:

(i) The trivial equilibrium point $E^*_0$ is globally stable if $r < \delta$ and $c_1 < \mu$.

(ii) The boundary equilibrium point $E_1^*$ is locally unstable if $\Gamma > 0$ and either $\dfrac{rmd}{\delta K}+c > \Gamma$ or $\dfrac{rmd}{\delta K}+c < \Gamma$.

(iii) The positive interior equilibrium $E^*_2 = \left(x^*, y^*, z^*\right)$ is locally stable if the Routh-Hurwitz stability criteria are satisfied.

Proof. The proof for the global stability for $E^*_0$ is similar to the proof of Theorem 2.3 (i) and is therefore omitted.

Evaluating $J^{**}$ at $E_1^*$ yields

A computation of the characteristic equation of the Jacobian $J_1^{**}$ is given by

Clearly, the real part of one of the eigenvalues $\lambda_1 = \dfrac{\alpha _1 y^*}{d_1+y^*}$ is positive and therefore the boundary equilibrium point $E_1^*$ is locally unstable.

We similarly evaluate the Jacobian $J^{**}$ at $E_2^*$ and obtain the following characteristic equation:

where

The positive interior equilibrium $E^*_2 = \left(x^*, y^*, z^*\right)$ is locally stable if $\theta_1 > 0$, $\theta_2 > 0$, $\theta_3 > 0$ and $\theta_1 \theta_2-\theta_3 > 0$ by the Routh-Hurwitz stability criteria.

3.   Bifurcation analysis

A bifurcation occurs when the behavior of a dynamical system changes when a parameter is varied. The said parameter at which this change occurs is known as the bifurcation parameter ^[28]. Bifurcation analysis is useful in measuring these qualitative changes and gives information on the parameters at which the system transitions from being stable to unstable and vice versa. We are interested in understanding the qualitative effects of the rate of cannibalism $c$ on the system.

3.1.   Co-dimension one bifurcations

3.1.1.   Hopf bifurcation

We shall present conditions under which the system (2.3) undergoes a Hopf bifurcation near the interior equilibrium point.

Theorem 3.1. If the rate of cannibalism $c$ in the juvenile prey population $x$ crosses the threshold value $c_H$, the system experiences a Hopf bifurcation around the interior equilibrium $E_2^*$ if the following conditions are satisfied:

and

Proof. We reconsider the characteristic Eq (2.13) of the form

for the occurrence of a Hopf bifurcation with roots $\lambda_{1}(c_{H}) = i \sqrt{\theta_{2}(c_{H})},$ $\lambda_{2}(c_{H}) = -i \sqrt{\theta_{2}(c_{H})},$ $\lambda_{3}(c_{H}) = - \theta_{1}(c_{H}) < 0$. Clearly, we have that $\theta_3(c_{H}) = \theta_1(c_{H})\theta_2(c_{H})$. We now need to establish the transversality condition

to verify the existence of periodic solutions bifurcating around $E_2^*$ at $c = c_{H}$. We substitute $\lambda_{k}(c) = \Gamma(c)+i\Lambda(c)$ into (3.3) and compute the derivative. We get

where

We solve for $\Gamma^{\prime}(c_{H})$ from the linear systems in (3.5) and (3.6) using Cramer's rule. We note that, at $c = c_{H},$ $\Gamma(c_{H}) = 0$ and $\Lambda(c_{H}) = \sqrt{\theta_{2}(c_{H})}$, resulting in

We now have

on condition that $\left[\theta_{1}(c)\theta_{2}(c)\right]^{\prime}_{c = c_{H}} - \theta_{3}^{\prime}(c_{H}) \neq 0.$

Therefore, the transversality condition is established, implying that the system experiences a Hopf bifurcation around $E_2^*$ at $c = c_{H}.$

Example 1. From Figure 1, we report that system (2.1) experiences a Hopf bifurcation around the point $E_2(x^*, y^*, z^*) = (0.1975, 0.39985, 0.366888)$ at $c = 0.442259$. With the help of Matcont software, the first Lyapunov coefficient is $\sigma = -6.66774046e^{-4}$ and hence the Hopf bifurcation is supercritical.

3.1.2.   Saddle-node bifurcation

The following theorem relates to the existence of a saddle-node bifurcation for the cannibalism rate $c$.

Theorem 3.2. The model (2.1) undergoes a saddle-node bifurcation around $E_2^*$ at $c = c^*$ when the conditions $det(J^{**}) = 0$ and $tr(J^{**}) < 0$ are satisfied by system parameters.

Proof. We apply the Sotomayor's theorem ^[27] in this proof to show the occurrence of a saddle-node bifurcation at $c = c^*$. At $c = c^*$, we can obtain $det(J^{**}) = 0$ and $tr(J^{**}) < 0$. This indicates that, $det(J^{**})$ admits a zero eigenvalue. Now let $U = (u_1, u_2, u_3)^T$ and $V = (v_1, v_2, v_3)^T$ be the eigenvectors of $J^{**}$ and $J^{**T}$ corresponding to the zero eigenvalue respectively.

We have that, $U = \left(\dfrac{-W}{B_1}, 1, \dfrac{-B_4}{B_5} \right)^T$ and $V = \left(\dfrac{-m}{A_1}, 1, \dfrac{-A_4}{A_5} \right)^T$ where $W = r\left(1-\dfrac{x^*}{K}\right)$, $A_1 = B_1 = \dfrac{-c x^* (2 d+x^*)}{(d+x^*)^2}+c_1-\mu -m-\dfrac{ry^*}{K},$ $B_4 = \dfrac{\alpha _1 d_1 z^*}{(d_1+y^*)^2},$ $A_4 = \dfrac{- \alpha y^*}{d_1+y^*}\; {\rm{and}}\; A_5 = B_5 = \dfrac{\alpha _1 y^*}{d_1+y^*}-2 \nu z^*$.

Furthermore, let $S = (S_1, S_2, S_3)^T$ where

Now,

and

Therefore, system (2.1) by Sotomayor's theorem undergoes a saddle-node bifurcation at $c = c^*$ around $E_2^*$.

Example 2. Again from Figure 1, a numerical experiment shows that system (2.1) undergoes two saddle-node bifurcations around $E_2(x^*, y^*, z^*) = (2.47282, 8.6393, 0.9260)$ at $c = 0.711034$ and around $E_2(x^*, y^*, z^*) = (0.45301, 0.8857, 0.5621)$ at $c = 0.21777$ respectively.

3.2.   Co-dimension two bifurcations

We further investigate the likely occurrence of co-dimension two bifurcations by performing a numerical continuation of co-dimension one bifurcations in system (2.1). We explore the $c-K$ space to understand the dynamics of the relationship between the rate of cannibalism attack $c$ and the juvenile prey carrying capacity $K$.

3.2.1.   Bogdanov-Takens bifurcation

Bogdanov-Takens bifurcation (BT) occurs when the critical equilibrium $(x^*, y^*, z^*)$ has a zero eigenvalue of multiplicity two after being evaluated at the Jacobian in Eq (2.11) in a two-parameter plane. We give a numerical example to illustrate a BT bifurcation in model (2.1) as observed in Figure 2. A BT bifurcation occurs for $(c, K) = (0.0564, 3.5405)$ around $E_2(x^*, y^*, z^*) = (0.6485, 1.3401, 0.6601)$. After the Jacobian in Eq (2.11) is evaluated at $E_2(x^*, y^*, z^*)$, the obtained eigenvalues are $\lambda_1 = -1.336089$, $\lambda_2 = \lambda_3 = 0$.

3.2.2.   Cusp bifurcation

A cusp bifurcation occurs when two curves at which points undergo a saddle-node bifurcation intersect tangentially. Cusp bifurcations lead to local bistability and hysteresis. This bistability result can be seen in Figure 1 for values of $c$ between $0.21777$ and $0.711034$ where there are two stable branches and an unstable branch of equilibrium points. We provide a numerical example of a cusp bifurcation in Figure 2. The cusp bifurcation is shown in the $c-K$ plane for system (2.1) around $E_2(x^*, y^*, z^*) = (0.8012, 1.7564, 0.7180)$ for $(c, K) = (0.03167, 3.29397)$.

4.   Numerical experiments

In this section, we provide numerical simulations to corroborate our theoretical findings. The parameters used in our simulations are auxiliary. We use MATHEMATICA $13.0$, MATCONT ^[29] and Python programming language in generating our plots. In Figure 3, numerical simulations show that the system is destabilized when cannibalism is absent. The system shows oscillatory dynamics and a limit cycle is formed for a set of chosen parameters. When cannibalism is introduced into the juvenile prey population at some rate, the oscillatory dynamics disappear and the system reaches a stable equilibrium as seen in Figure 4. The presence of cannibalism therefore has a stabilizing effect. Simulations in Figure 5 show how cannibalism can stabilize and destabilize the system. When there is no cannibalism, the system is stable at the coexistence equilibrium $E_2(x^*, y^*, z^*) = (6.28, 27.3234, 0.97537)$. As the rate of cannibalism increases to some $c^* = 0.442$, the system becomes destabilized and a stable limit cycle is formed. A further increase in the cannibalism rate brings the system back to a stabilized state at $E_2(x^*, y^*, z^*) = (0.16013, 0.33479, 0.32669)$ when $c^* = 0.6$. Thus we observe that an increase in the rate of cannibalism leads to a decrease in the population densities of the prey and predator.

We also validate our findings on the global stability of the trivial steady state of the system when cannibalism is absent in Figure 6(i) when the juvenile prey growth rate $r$ is less than the death rate $\delta$ of the adult prey. In Figure 6(ii), we also show that when cannibalism is present and if both restrictions $r < \delta$ and $c_1 < \mu$ hold, the trivial steady is also globally stable. We see another interesting dynamics under a parametric regime that cannibalism at a certain rate can induce coexistence while the absence of it leads to extinction of all the populations in Figure 7.

5.   Subsystem in the absence of predators

We want to study the dynamics of system (2.1) when there are no predators. Our interest is to explore if the presence of predators has an effect in the generation of cyclic dynamical patterns in the system. We shall only compute the equilibria in the Appendix 6 and refrain from performing stability analysis on the subsystem. Thus system (2.1) reduces to

We next show that in the absence of predators, the subsystem (5.1) cannot produce oscillatory dynamics via the following theorem:

Theorem 5.1. The system (5.1) has no limit cycles for $c > c_1 > 0$.

Proof. We use the Dulac theorem to show that limit cycles do not exist in system (5.1). We consider the following Dulac function

where both $x$ and $y$ are non-zero. Then, we have

Hence there is non-existence of limit cycles in system (5.1).

Corollary 1. The system (5.1) has no limit cycles for $c = c_1 = 0$.

6.   Discussion and conclusions

In this current manuscript, a three system ODE stage-structured predator-prey model is considered with cannibalism occurring only in the juvenile prey. Cannibalism is a trophic interaction and its impact can alter the relationship between predators and prey as it is known to occur in over 1300 species ^[1]. Our motivation to study the effect of cannibalism stems from the report on cannibalism occurring among tadpoles in a specie of treefrogs ^[21] and also from the study by Aladeen et al. ^[13] that there are no mathematical models that explore cannibalism in juvenile prey.

In the absence of cannibalism, our simulations show that the prey and predator populations exhibit oscillatory dynamics for a chosen set of parameters as seen in Figure 3. An increase in the cannibalism attack rate changes the coexistence point from an unstable to a stable state as seen in Figure 4. This shows that, cannibalism has a stabilization effect based on our parameter values and supports the findings of ^[30]. Also under a certain parameter regime, we report that cannibalism can have both a stabilization and destabilization effect as the rate of attack is increased/decreased and corroborates findings in ^[11]. This leads to a decrease in the population density of both prey and predator populations as seen in Figure 5. The decrease in the steady state of the populations as a result of cannibalism is an observed phenomena, for example, in cyclopoid copepod simulation ^[31]. This result pertaining to population size decrease also has applications to invasive species control. It has been reported in the Florida Everglades that prey populations are severely declining as a result of the presence of invasive Burmese pythons ^[32]. Results in Figure 5 give insight on cannibalism being able to help maintain prey populations at lower levels, causing a decline in the invasive predator population through intraspecific competition. In this case, we see that juvenile prey cannibalism plays the role of population regulation.

Also, numerical results seen in Figures 6 and 7 show that, in the absence of cannibalism, all populations go extinct. Interestingly, the presence of cannibalism can also cause all populations to go extinct or coexist for a parametric regime. It is this coexistence which is labelled the so called "life boat mechanism". A biological implication of this coexistence is that survivors of cannibalism are better fed and are competent in giving rise to new offspring when they mature into adult prey ^[1]. We also observed that the stage-structured system possessed rich bifurcation dynamics. The system was found to experience Hopf, saddle-node, Bogdanov-Takens and cusp bifurcations respectively for various parameter regimes via numerical simulations. From an ecological point of view, the occurrence of a Hopf bifurcation indicates that in a stable predator-prey environment, cannibalism in the juvenile prey can cause oscillations in the population and drive the instability of the system. The numerical example of the saddle-node bifurcation indicates that the system (2.1) exhibits bistability for a chosen set of parameters and hence sensitive to initial conditions. In this case, the populations will continue to coexist. We leave the rigorous proofs for the occurrence of cusp and Bogdanov-Takens bifurcations as a future work. See ^[33] for more complex two parameter bifurcations. In all, these results shed light on the role juvenile prey cannibalism plays in the dynamics of the populations.

Furthermore, we established that the subsystem model cannot exhibit oscillatory dynamics via Dulac's criterion. Thus the presence of predators and the Holling type Ⅱ functional response can be the catalyst for the observance of cycling dynamics in the three system ODE stage-structured model. It will be worth exploring different functional responses in modeling cannibalism terms and adult prey and predator interactions for interesting dynamics in future. It is still interesting to note that, there is no consensus to the impacts of cannibalism on dynamical systems with regards to its stabilization and destabilization effects. Therefore, future work in both mathematics and biology is needed to fully understand the complex consequences that cannibalism plays in population dynamics. Future work herein may also include incorporating a time delay in the maturation of juveniles to adults in the prey population and study the overall system dynamics. Fractional predator-prey models have been studied recently ^[34]. The dynamics of fractional stage-structured model with cannibalism remains unexplored. It will also be interesting to further extend the ODE model to a spatial model and explore the possible occurrence of Turing patterns in one and two dimensions which has applications in biocontrol.

Acknowledgements

EMT, KC, AD and CM would like to acknowledge valuable support from the National Science Foundation via grant number 1851948.

Conflict of interest

The author declare that they have no conflict of interest.

Appendix

A.1

The coefficients seen in Eq (2.10) are given as:

$\eta_5 = -\dfrac{\delta \nu \left(\nu \left(\delta K (c+\mu -c_1)+d m r+4 \delta d_1 r+K m (\delta -r)\right)+2 \alpha \alpha _1 r\right)}{K}$,

$\eta_4 = \dfrac{1}{K}\left[H_1+ \nu ^2 \left(H_2+H_3-6 \delta ^2 d_1^2 r\right)-\alpha ^2 \alpha _1^2 r \right]$,

where

$\eta_3 = \dfrac{1}{K}\left[H_4+2 d_1 \nu ^2 \left(H_5+H_6\right)-\alpha ^2 \alpha _1^2 K \left(c-c_1+\mu +m\right)\right]$,

where

$\eta_2 = -\dfrac{1}{K}\left[H_7+d_1 \nu \left(4 \delta d_1 \left(\delta K (c+\mu -c_1)+d m r+K m (\delta -r)\right)+H_8\right)\right]$,

where

$\eta_1 = -\dfrac{d_1^2 \nu \left(H_9+4 d d_1 K m \nu \left(-c_1 \delta +\delta \mu +\delta m-m r\right)\right)}{K}$,

where

$\eta_0 = d d_1^4 m \nu ^2 \left(c_1 \delta -\delta \mu +m (r-\delta)\right)$.

A.2. Equilibria when cannibalism is absent

In this subsection, we find that the subsystem (5.1) has two equilibria:

${\rm{(a)}} \ P_0 = (0, 0)$ and

${\rm{(b)}}$ a unique coexistence equilibrium $P_1 = (x^*, y^*)$, where

$P_1 = (x^*, y^*)$ exists if $(r-\delta) > \dfrac{\delta \mu}{m} > 0$ provided that $\delta < r$.

A.3. Equilibria when cannibalism is present $(c > c_1 > 0)$.

By solving $\dfrac{dx}{dt} = \dfrac{dy}{dt} = 0$ in subsystem (5.1), we obtain a trivial equilibrium and a unique interior equilibrium given by

$({\rm{i}})\ P^*_0 = (0, 0)$,

$({\rm{ii}})$ $P^*_1 = (x_1^*, y_1^*),$

$({\rm{ii}})$ $P^*_2 = (x_2^*, y_2^*),$ where

where

We note that both $P^*_1$ and $P^*_2$ exist if $x_1^* > 0$ and $x^*_2 > 0$.