Sex-biased predation and predator intraspecific competition effects in a prey mating system

1.   Introduction

Ecologists and applied mathematicians have extensively investigated the role of predation by predators in controlling the numbers of their prey ^[1,2,3,4,5]. One type of predation that has been observed is sex-biased predation. Sex-biased predation is the degree by which both sexes of prey are predated disproportionately to their abundance or general ratio ^[6]. Often, male and female prey differ in vulnerability. This can be due to differing factors such as size, speed, and age, often referred to as sexual dimorphism. Sexual dimorphism can be in part due to the different responsibilities associated with reproduction and childbearing in females and the role of males in different species of prey ^[7]. Many species assign different roles to male and female prey populations, such as foraging for food, nursing their young, and staking territory ^[6]. These roles may result in one sex becoming more exposed to their predators than the other. There are also predation risks associated with mate-finding in both sexes of prey depending on who searches and who signals ^[8]. It is well known that in aquatic environments invasive species cause extensive damage to wetland biodiversity and ecosystem function ^[9,10,11]. For example, the freshwater snail Pomacea Canaliculata is known to be one of the world's worst invasive species by the Global Invasive Species Programme ^[12]. They grow very fast, have high fecundity and cause extensive damage to agricultural production and wetland biodiversity ^[7]. Xu et al. investigated and reported the existence of male-biased predation in Pomacea Canaliculata when introduced to turtles Chinemys Reevesii as predators. The male-biased predation could be attributed to males traveling longer distances in search for mates over larger areas when the density of snails is low. Hence they have a higher encounter rate with predators. Other prey species experiencing male-biased predation include guppy fish ^[13] and male possums ^[14]. Female-biased predation has also been observed in mosquito fish Gambusia affins ^[15]. Evidence from the fieldwork of Britton and Moser ^[16] showed that captive herons, predators of Gambusia, almost exclusively preyed on females even when heavily outnumbered by males.

Sex-biased predation can have effects on the life history traits of species. Research experiments conducted by DeGabriel et al. ^[14] revealed that about $60\%$ of marked common brushtail Trichosurus vulpecula male possums in tropical populations were preyed on by pythons due to their mobility when searching for mates. Their results also brought to light that male possums showed traits of investing in increased growth during their early stages of life so as to maximize their potential of mating. The introduction of some prey species to predation threat can result in behavioral adjustments such as changing reproductive strategies or courting behavior. For instance, male guppies Poecilia reticulata are more favored to be involved in sneak copulations rather than courtship when exposed to predators ^[17]. This leads to missed mating opportunities and negatively affects male reproductive success ^[18,19]. In all, sex-biased predation affects the sex ratio and the evolution of a population. The consequences of an imbalanced sex-ratio include competition for mates ^[20], parental care ^[21], mate choice ^[22], alternative mating strategies ^[23] etc. Boukal et al. ^[24] studied a predator-prey system where predators profited from prey exhibiting sexual dimorphism traits. Their findings showed that sex-biased predation depended on the interaction between the prey mating dynamics and predation bias. Also, male or female biased predation and mate-finding Allee effect had an impact in (de)stabilizing the prey mating system. However, there is still more to be explored on the impact of sex-biased predation on population stability and persistence ^[24].

In many ecological communities, competition is a common interaction observed when individuals vie for scarce resources and habitat space. When the density of a population is high, intraspecific competition negatively affects survivorship and/or fecundity ^[25]. The impacts of predator intraspecific competiton on invasive pests have been less understood ^[26,27,28]. However, research done by Parshad et al. ^[4] showed that competition is able to prevent unbounded growth of an introduced predator to drive a target pest population to extinction when supplied with high quality and sufficient additional food resource. Furthermore, intraspecific competition effects enhanced biocontrol and also led to the occurrence of non-standard co-dimension one and two bifurcation structures ^[4]. In this work, we shall incorporate predator intraspecific competition to our modeling framework. As mentioned earlier, the impact of sex-biased predation in mating models has been less investigated ^[24]. The aim of this paper is to understand the long term dynamics and consequences of sex-biased predation and the effects of predator intraspecific competition from a mathematical modeling standpoint by investigating the following ecological questions:

● How does sex-biased predation affect the sex-ratio of the prey?

● Can sex-biased predation stabilize or destabilize the population dynamics?

● How does predator competition affect the overall population dynamics?

● What are the ecological implications of the bifurcation structures obtained in the study?

● Does predator presence affect the population dynamics?

We organize the rest of the paper as follows: In Section 2 we formulate our mathematical model with its underlying ecological assumptions. We provide preliminary results such as positivity and boundedness of solutions to our model in Section 3. We study a special case of our system when predators are not introduced in Section 4. We show that the system cannot produce oscillatory dynamics when predators are absent. In Section 5 we provide analytical guidelines on equilibrium points and their stability properties when the predator population is present. We discuss local bifurcation results in Section 6 and observe that sex-biased predation can have destabilizing and stabilizing effects. We provide numerical simulations to validate our theoretical findings in Section 7. We notice that the sex-ratio can be constant and imbalanced, constant and balanced or cyclic which has ecological consequences. See Figures 4$(c)$, 6$(c)$ and 5$(c)$, respectively. We give a discussion on our findings in Section 8 and a conclusion in Section 9.

2.   Model formulation

Here, we consider a sex-structured prey population divided into female and male classes, who are depredated based on sex-bias by an unstructured predator population. We let $f, m$ and $p$ denote the female, male and predator population, respectively, at any time instant $t$. We make the following assumptions in the formulation of our model:

(ⅰ) The total prey population can grow up to a carrying capacity $K$, which is determined by limited resources and so we include a logistic term $L$ to modulate growth. Also, the prey population growth is proportional to interactions between the males $m$ and females $f$ ^[29].

(ⅱ) We assume that there is no bias in the sex-ratio of the prey at birth.

(ⅲ) Predation is sex-biased. We let $r$ denote the sex-biased predation rate on the male population. As a result, $1-r$ is the sex-biased predation rate on females. Also, $0 < r < 1$.

(ⅳ) We use the Holling type Ⅱ functional response to describe the relationship between the predator and its prey. This is described by the $\Phi_1(f, m)$ and $\Phi_2(f, m)$ terms, where $\alpha$ is the attack rate of the predator and $\nu$ is the time the predator uses in handling its prey.

(ⅴ) We assume natural death rates $\delta$ for the prey and $\delta_{1}$ for the predator.

(ⅵ) We also assume there is intraspecific competition between the predators for resources, which is modeled by the $-cp^2$ term.

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

where

with positive initial conditions $f(0) = f_0$, $m(0) = m_0$ and $p(0) = p_0$. L is the logistic term. All parameters used are assumed positive and their descriptions are given in Table 1.

3.   Preliminary results

This section provides basic results on the positivity and boundedness of solutions for system (2.1).

3.1.   Positivity of solutions

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

Lemma 3.1. Consider the following system of ODEs:

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

if and only if

we have

3.2.   Boundedness of solutions

We state the following theorem:

Theorem 3.2. All solutions of system (2.1) which initiate from $\mathbb{R}_3^+$ are bounded.

Proof. The female and male populations are already assumed to grow up to a carrying capacity $K$. So we have $f\leq K$ and $m\leq K$. Therefore

by the comparison principle. Next, simple calculations show that ${\rm{lim sup}}\; p(t) \leq \dfrac{\alpha K(\gamma_1+\gamma_2r)}{c}$ as $t \rightarrow \infty$. Hence, all solutions which initiate from $\mathbb{R}_3^+$ are bounded for system (2.1). □

4.   Case where there are no predators

In the absence of predators, system (2.1) reduces to a classic mating system given by

where $L = 1-\dfrac{f+m}{K}$ with positive initial conditions $f(0) = f_0$, $m(0) = m_0$. Results pertaining to the equilibrium points of system (4.1) and their stability are well studied ^[32,33]. We will prove the non-existence of limit cycles for system (4.1) using the Dulac theorem.

Theorem 4.1. The system (4.1) cannot exhibit limit cycle dynamics.

Proof. We apply the Dulac theorem to show that system (4.1) does not exhibit cyclic dynamics. Consider the Dulac function

where both $f$ and $m$ are non-zero. Then, we have

Hence, system (4.1) cannot produce limit cycle dynamics. □

5.   Case where predators are present

5.1.   Equilibria

To obtain the equilibria for system (2.1), we solve $F(f, m, p) = 0$, $M(f, m, p) = 0$ and $P(f, m, p) = 0$ simultaneously. The system possesses the following non-negative equilibria:

(a) $E_{0} = (0, 0, 0)$,

(b) $E_{1} = (f_{1}^{*}, m_{1}^{*}, 0)$ where $f_1^* = m_1^* = \dfrac{K}{4} (1\pm \sqrt{1-\Gamma})$ and $\Gamma = \dfrac{16\delta}{K\beta}.$ We note that when $\Gamma = 1$, $f_{1_\pm}^* = m_{1_\pm}^* = \dfrac{K}{4}$. When $\Gamma > 1$, there are no real equilibrium points, and when $\Gamma < 1$, there are two distinct interior equilibria given by $f_{1_+}^* = m_{1_+}^* = \dfrac{K}{4} (1+ \sqrt{1-\Gamma})$ and $f_{1_-}^* = m_{1_-}^* = \dfrac{K}{4} (1- \sqrt{1-\Gamma})$.

(c) $E_2 = (f_2^*, m_2^*, p_2^*)$ where $p_2^* = \dfrac{1}{c}\left[\dfrac{\gamma_1 \left(1-r\right)\alpha f_2^*+\gamma_2 r \alpha m_2^*}{1+\nu\left(f_2^* +m_2^* \right)}-\delta_1 \right]$.

Now, substituting $p_2^*$ into $F(f, m, p) = 0$ yields

where $B_1 = \gamma_1 \left(1-r \right)\alpha f_2^*+\gamma_2 r \alpha m_2^*-\delta_1 \left[1+\nu \left(f_2^*+m_2^* \right) \right]$.

The real positive roots $f_2^*$ for the solution to equation (5.1) will be in terms of $m_2^*$ and are substituted into $M(f, m, p) = 0$ to find an explicit expression for $m_2^*$. However, it is difficult to find this expression due to very complicated algebraic calculations. We will illustrate the existence and stability properties of $E_2$ via time series simulations and bifurcations when certain control parameters are varied.

5.2.   Local stability analysis

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

where

Theorem 5.1. The extinction state $E_{0}$ is locally stable.

Proof. We evaluate $J^*$ at $E_0$ and obtain

The eigenvalues are $\lambda_1 = -\delta < 0, \; \lambda_2 = -\delta < 0, \; \lambda_3 = -\delta_1 < 0.$ Therefore the extinction state $E_{0}$ is locally stable. □

Similarly, we evaluate $J^*$ at $E_1$ and obtain

where $\rho = \alpha \gamma _1 f^*(1-r)+\alpha \gamma _2 m^* r$.

Theorem 5.2. If $\Gamma = 1$, then $E_{1} = (\frac{K}{4}, \frac{K}{4}, 0)$ is locally stable provided $\beta < \dfrac{16 \delta}{K}$ and $\dfrac{\alpha K \left[ \gamma_1 (1-r) + \gamma_2 r \right]}{2(K \nu +2)} < \delta_1.$

Proof. Let $\Gamma = 1$, $\beta < \dfrac{16 \delta}{K}$ and $\dfrac{\alpha K \left[ \gamma_1 (1-r) + \gamma_2 r \right]}{2(K \nu +2)} < \delta_1.$ Then, evaluating $J^*$ at $E_1$ yields

The characteristic equation for $J^*_{{E_1}_{(\Gamma = 1)}}$ is given by

and the eigenvalues are $\lambda_1 = \dfrac{K \beta}{16}-\delta, \; \lambda_2 = -\delta < 0, \; \lambda_3 = \dfrac{-2 \delta _1 (K \nu +2)+\alpha \gamma _1 K(1-r)+\alpha \gamma _2 K r}{2(K \nu +2)}.$ Since $\Gamma = 1$, $\beta < \dfrac{16 \delta}{K}$ and $\dfrac{\alpha K \left[ \gamma_1 (1-r) + \gamma_2 r \right]}{2(K \nu +2)} < \delta_1$ by our assumption, $E_1$ is locally stable. □

Now, when $\Gamma < 1$, we have the characteristic equation

where $A_1 = -\eta_1-\eta_4-\eta_5$, $A_2 = \eta_4 \eta_5+\eta_1(\eta_4+\eta_5)-\eta_2 \eta_3$ and $A_3 = \eta_2 \eta_3 \eta_5-\eta_1\eta_4 \eta_5$ with $\eta_1 = \dfrac{\beta m^* (K-2 f^*-m^*)}{2 K}-\delta$, $\eta_2 = \dfrac{\beta f^* (K-f^*-2 m^*)}{2 K}$, $\eta_3 = \dfrac{\beta m^* (K-2 f^*-m^*)}{2 K}$, $\eta_4 = -\delta +\dfrac{\beta f^* (K-f^*-2 m^*)}{2 K}$ and $\eta_5 = \dfrac{-\delta _1 (\nu (f^*+m^*)+1)+\rho}{\nu (f^*+m^*)+1}$.

The equilibrium point $E_1 = (f_1^*, m_1^*, 0)$ is locally stable if $A_1 > 0, \; A_2 > 0, \; A_3 > 0, \; {\rm{and }}\; A_1 A_2- A_3 > 0$ by the Routh-Hurwitz theorem.

In ascertaining the local stability of the coexistence equilibrium $E_2$, the characteristic equation of $J^*_{E_2}$ is

where

By applying the Routh-Hurwitz stability criteria, $E_2$ is asymptotically stable if $\sigma_1 > 0, \; \sigma_2 > 0, \; \sigma_3 > 0, \; {\rm{and }}\; \sigma_1 \sigma_2- \sigma_3 > 0$.

5.3.   Global stability analysis

Theorem 5.3. The extinction state $E_0$ is globally stable if $\beta < \dfrac{2 \gamma_1 \delta}{K(\gamma_1+\gamma_2)}.$

Proof. Suppose that $\beta < \dfrac{2 \gamma_1 \delta}{K(\gamma_1+\gamma_2)}$ and consider the Lyapunov function $V(t) = \gamma_1 f(t)+\gamma_2 m(t)+p(t)$ where $\gamma_1$ and $\gamma_2$ are positive as already assumed. Clearly, $V = 0$ at $(f, m, p) = (0, 0, 0)$ and $V > 0$ when $(f, m, p) \neq (0, 0, 0).$ Now, we compute the time derivative of V and get

Therefore $E_0$ is globally stable. Hence the proof is complete. □

6.   Bifurcation analysis

Bifurcation studies provide insights on the qualitative changes of the behavior of a dynamical system when one or more parameters are varied. When a bifurcation happens, the stability attributes of equilibrium points and periodic orbits changes. Of particular interest to us is the sex-biased predation rate parameter $r$, the intraspecific competition among predators parameter $c$ and the predator attacking rate $\alpha$.

6.1.   Hopf bifurcation

A Hopf bifurcation is said to occur when there is an appearance or disappearance of a periodic orbit when there is a change in the local stability of an equilibrium point. The following theorem relates to the occurrence of a Hopf bifurcation for the sex-biased predation rate $r$.

Theorem 6.1. If the sex-biased predation rate $r$ crosses a threshold value at $r = r_H$, system (2.1) experiences a Hopf bifurcation around the coexistence state $E_2$ if the following conditions hold:

and

Proof. Let us suppose that the characteristic equation (5.3) is of the form

with roots $\lambda_{1}(r_{H}) = i \sqrt{\sigma_{2}(r_{H})},$ $\lambda_{2}(r_{H}) = -i \sqrt{\sigma_{2}(r_{H})},$ $\lambda_{3}(r_{H}) = - \sigma_{1}(r_{H}) < 0$. Clearly, $\sigma_3(r_{H}) = \sigma_1(r_{H})\sigma_2(r_{H})$. The next step is to validate the transversality condition

to show that periodic solutions exist and bifurcates around $E_2^*$ at $r = r_{H}$. We substitute $\lambda_{j}(r) = \xi (r)+i\Lambda(r)$ into (6.3) and compute the derivative. We obtain

where

We apply the Cramer's rule to solve for $\xi^{\prime}(r_{H})$ in the linear systems in (6.5) and (6.6). Observe that at $r = r_{H},$ $\xi(r_{H}) = 0$ and $\Lambda(r_{H}) = \sqrt{\sigma_{2}(r_{H})}$, which yields

Simple calculations show that

subject to $\left[\sigma_{1}(r)\sigma_{2}(r)\right]^{\prime}_{r = r_{H}} - \sigma_{3}^{\prime}(r_{H}) \neq 0.$

This establishes the transversality condition and hence the occurrence of a Hopf bifurcation around $E_2$ at $r = r_{H}.$ □

6.2.   Saddle-node bifurcation

A saddle-node bifurcation occurs when two equilibria collide and annihilate each other. We shall use Sotomayor's theorem ^[34] to show that system (2.1) experiences a saddle-node bifurcation at a critical intraspecific competition threshold value $c = c^*$.

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

Proof. We use Sotomayor's theorem ^[34] to show that system (2.1) experiences a saddle-node bifurcation at $c = c^*$. At $c = c^*$, we can have $det(J^{*}) = 0$ and $tr(J^{*}) < 0$. This shows that $J^{*}$ admits a zero eigenvalue. Define $X = (x_1, x_2, x_3)^T$ and $Y = (y_1, y_2, y_3)^T$ to be the nonzero eigenvectors of $J^{*}$ and $J^{*T}$ corresponding to the zero eigenvalue, respectively. For $J^*_{23} = \dfrac{J^*_{22}J^*_{13}}{J^*_{12}}$ and $J^*_{33} = \dfrac{J^*_{32}J^*_{13}}{J^*_{12}}$ with $J^*_{12}\neq 0$ and $J^{*T}_{32} = \dfrac{J^{*T}_{22}J^{*T}_{31}}{J^{*T}_{21}}$ and $J^{*T}_{33} = \dfrac{J^{*T}_{23}J^{*T}_{31}}{J^{*T}_{21}}$ with $J^{*T}_{21}\neq 0$, we have

$X = \left(0, -\dfrac{J^*_{13}}{J^*_{12}}, 1 \right)^T$ and $Y = \left(0, -\dfrac{J^{*T}_{31}}{J^{*T}_{21}}, 1 \right)^T$.

Furthermore, let $Z = (Z_1, Z_2, Z_3)^T$ where

Now,

and

Therefore, system (2.1) by Sotomayor's theorem experiences a saddle-node bifurcation at $c = c^*$ around $E_2$ and the proof is complete. □

Theorem 6.3. If the predator attack rate $\alpha$ crosses a threshold value at $\alpha = \alpha_H$, system (2.1) experiences a Hopf bifurcation around the coexistence state $E_2$ if the following conditions hold:

and

Proof. The proof is similar to Theorem 6.1 and is therefore omitted. □

7.   Numerical experiments

In this section, we provide numerical simulations to support our theoretical findings. We used the Python programming language to generate our time series and phase plots. Figure 4 shows the existence of three biologically feasible equilibria for system (2.1) for a chosen set of parameter values. Figure 4$(a)$ shows the extinction state of all the populations, $(b)$ shows the predator free state with equilibrium point $E_1 = (19.9328, 19.9328, 0)$ and $(c)$ shows the coexistence state of all populations with equilibrium point $E_2 = (10.1991, 23.2547, 5.3464)$. It is worth noting that the sex-ratio is constant and balanced from Figure 4$(b)$, and constant but imbalanced in 4$(c)$. In Figure 5, numerical simulations show oscillatory dynamics for a different parameter set. Similarly, we observe that the sex-ratio is not constant but cyclic as seen in Figure 5$(c)$. However, the sex-ratio is constant and balanced when there is no sex-biased predation and the population is cyclic as seen in Figure 6. We validated a sufficient condition in Theorem 5.3 with regards to the global stability of the extinction state with experiments seen in Figure 7. Thus when the birth rate $\beta$ is less than $\frac{2 \gamma_1 \delta}{K(\gamma_1 + \gamma_2)}$, all species die out irrespective of the initial population density.

For generating our bifurcation plots, we used MATLAB version R2019a and MATCONT ^[35] software. Figure 1 shows the existence of multiple Hopf bifurcations for the sex-biased predation parameter $r$ when all other parameters are fixed. When $r$ is increased, the coexistence equilibrium loses its stability at the critical threshold value $r^* = 0.4098569$ and oscillatory dynamics emerge. At this threshold value, there is an occurrence of a Hopf bifurcation around the coexistence equilibrium $E_2 = (10.146745, 14.53349, 6.92543)$. The first Lyapunov coefficient is computed with MATCONT and is given by $\chi = 2.90744853e^{2}$ and hence a subcritical Hopf bifurcation. A further increase in $r$ leads to the occurrence of another Hopf bifurcation at the threshold value $r^* = 0.590143$ around $E_2 = (14.53348, 10.14674, 6.92543)$. A similar computation of the Lyapunov coefficient gives $\chi = 2.90742685e^2$ and hence the Hopf bifurcation is subcritical. Once $r$ is increased again from $r^* = 0.590143$, the system regains its stability. This demonstrates that sex-biased predation has both stabilizing and destabilizing effects under certain parametric choices.

Numerical experiments also show that the system experiences a Hopf bifurcation for the predator attack rate $\alpha$ at $\alpha^* = 0.27503$ around $E_2 = (8.37687, 12.49559, 6.49191)$ as shown in Figure 2 with Lyapunov coefficient $\chi_1 = -1.925015e^{-4}$. The Hopf bifurcation is supercritical. Figure 3 also shows the occurrence of a saddle-node bifurcation for the intraspecific competition parameter $c$ at the critical value $c^* = 0.18428$ around $E_2 = (16.72081, 4.28190, 10.57134)$ where two coexistence equilibria collide and annihilate each other.

8.   Discussion

Our results show that sex-biased predation can cause an imbalanced sex-ratio in the prey population. See Figure 5(c) and compare with Figure 6(c). This has ecological consequences such as pairing behavior; male-male, male-female, female-female ^[36] and competition for mating access, which has implications for the success of natural populations. When the ratio of females to males is low, it also impacts population growth and hence has consequences for population dynamics, risk of extinctions and biodiversity conservation ^{[37,38,39,40]}. In the absence of predators, system (2.1) reduces to the classic mating system. We rigorously proved that the classic mating system cannot exhibit oscillatory dynamics via Theorem 4.1. However, the introduction of predators and the choice of a Holling type Ⅱ functional response can cause the populations to fluctuate via oscillations when the predation is sex-biased. See Figure 5. It is also possible to see cyclic dynamics even when there is no sex-biased predation ($r = 0.5$) as shown in Figure 6 as the population densities of males and females are equal and the sex-ratio is 1. Therefore predator presence affects the dynamics of the populations.

Under certain parametric regimes, we report that sex-biased predation can have stabilizing and destabilizing effects as the rate of the biased predation is increased or decreased. This corroborates findings in ^[24]. This is seen in Figure 1 when the stable interior equilibrium point loses its stability via a Hopf bifurcation and subsequently gains stability via another Hopf bifurcation as the sex-biased predation parameter $r$ is varied. Similar bifurcation results are seen in Figure 2 for the rate $\alpha$ at which predators attack prey. This shows that the populations can continue to thrive or persist. It is clear that increasing $r$ leads to an increase in the female population density and a decrease in the male population density as seen in Figure 1 and vice-versa. Depending on the severity of the biased predation, it can have implications for species conservation and extinction risks. Dynamically, when the predator intraspecific competition rate crosses some threshold $c = c^*$, a saddle-node bifurcation occurs when two coexistence states collide and annihilate each other, see Figure 3, and leaves the extinction state to be globally attracting. Our saddle-node bifurcation result is very useful in the context of biological control for those invasive species that are particularly susceptible to sex-biased predation. Examples include freshwater snail ^[12], mosquito fish ^[15] and guppy fish ^[13]. In 2017, the financial cost of controlling and mitigating the impact of invasive species in Australia was estimated to be ＄298.58 billion ^[41]. Also, the cost incurred in invasive species control in the United States was reported to total ＄1.2 trillion in 2020 ^[42]. It is therefore a paramount concern and curbing their damages to ecosystems and the shattering of fragile food webs will be of economic and societal value. Therefore, the introduction of natural enemies (predators) as agents of biocontrol can help in their eradication. Other mathematical approaches to invasive species control can be seen in ^[43]. From our system, this application is possible when competition for resources (prey) is less among predators as this has less negative impact on their growth/reproduction rate and their density in combating the invasive pest. Hence, there is a greater chance of eliminating the invasive pests. Conversely, the invasive pest will not go extinct when rate of intraspecific competition is high.

9.   Conclusions

In conclusion, our study showed that sex-biased predation has an impact on the dynamics of populations with regards to species coexistence. It can cause the population to stabilize as seen in Figure 4 and also fluctuate for a chosen parameter set as illustrated in Figure 5 showing a periodic solution and the periodic nature of the male to female sex-ratio. This has biological consequences for population dynamics ^{[20,21,22,23]} as already highlighted in the introduction section. The study also showed that the presence of predators can lead to cyclic population dynamics. Rich dynamical structures were also revealed via Hopf and saddle-node bifurcations and global stability results were obtained for the extinction state with application for invasive species control. See Figure 7 and proof of Theorem 5.3. In all, further investigation is needed to gain more insight on the impact of sex-biased predation and sexual dimorphism in prey on predator-prey interactions, as these could lead to very interesting and richer dynamics. It will be an interesting future work to consider other well known functional responses as well as mating functions in the modeling framework of system (2.1) and approaches from recent developments in mathematical modeling such as fractional-order density dependent models ^[44]. Finally, this work adds to the growing literature in investigating and considering the impacts of sex-biased predation to get a better understanding of evolutionary ecology dynamics ^[45].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

EMT, CO, MC and NK would like to acknowledge valuable support from the National Science Foundation via grant number 1851948 at Ursinus College.

Conflict of interest

The authors declare that they have no conflict of interest.