Qualitative analysis and new traveling wave solutions for the stochastic Biswas-Milovic equation

1.   Introduction

Classical competition theory predicts that two competing species (or even sub-populations of the same species) whose niches overlap perfectly cannot coexist, and one must competitively exclude the other ^[1]. The classical Lotka-Volterra ODE competition model (LVM) predicts this competitive exclusion dynamics, as well as co-existence and bi-stability, finding numerous applications ^[2,3,4,5,6]. Classical model literature models the competition by considering growth/fitness rates as inter and intra-species competition ^[2,7,8]. The growth or "fitness" term in the model is essentially a measure of reproductive success ^[9]. To this end, one can define "relative fitness" as a measure of differential reproductive success—that is, the proportion of the next generation's gene pool coming from a particular organism in contrast to other competing organisms ^[9]. All other parameters being equal in LVM, a difference in growth rates would point to differential fitness among the two competitors. A well-studied area within competition theory is that of adaptive dynamics ^[10]. Herein, the basic concept is that evolution occurs on a slower timescale than ecology. Once in a while, a mutation arises, and a small new population of mutants develops, which is identical to the existing population except for some specific trait. If that population can invade the existing population, then the new trait can persist and may take over the existing population ^[11]. A trait is evolutionarily stable if no small population with a different trait can invade it ^[12]. An important question in adaptive dynamics is: Can a small but efficient weaker sub-population invade a stronger sub-population? We aim to address this question in part in this manuscript. Invasive species are most successful in environments lacking close relatives ^[13,14]. This has been confirmed in experiments, where invasion success in microbial communities increases as phylogeny or "relatedness" between invader and invadee decreases ^[15]. Thus, an invasion is most likely to be successful if the invasive species faces low intraspecies competition whilst being a superior interspecies competitor ^[16,17]. However, in reality, a successful invasion results from the interaction of a myriad of factors, biological, environmental, landscape, and temporal - this is quantified via the concept of a species ecological niche ^[8,18]. A niche space can be thought of as a continuum, similar to $\mathbb{R}^{4}_{+}$, where the four axes are physical and environmental factors (say climate), biological factors (say resources), space, and time. A species niche is the response it has to each point in the space and the effect it has at each point ^[8,19]. The niche-based hypothesis posits invasive species are dexterous at using unexploited resources - that is, filling vacant niches or broadening their niche breadth if the opportunity presents itself ^[19,20]. One such "opportunity" is presented by climate change. Climate change is potentially as large a threat to biodiversity as habitat degradation ^[21]. Researchers focus on range shifts following local colonization and/or extinction of species. Herein, the "niche" concept resurfaces, in that species typically survive in a certain window of environmental/climatic conditions, defining the "niche space" so that if the environment changes, species will have to adapt, disperse, or otherwise be driven to extinction ^[22]. This becomes even more intriguing from an invasive species standpoint when climatic conditions turn favorable for an invader (possibly mutant or even a weaker sub-population) and unfavorable for the resident (or stronger sub-population) ^[23]. Here, we focus in part on a particular niche change pertaining to the soybean aphid Aphis glycines (Hemiptera: Aphididae).

The soybean aphid, Aphis glycines (Hemiptera: Aphididae) is an invasive pestiferous species, native to Asia, and was first discovered in the US in 2000 in Wisconsin ^[24]. It has a heteroecious holocyclic life cycle. In the spring, aphids emerge and produce three asexual generations on common buckthorn, Rhamnus cathartica L., then migrate to soybeans Glycine max L. Aphids continue to reproduce asexually on soybeans, producing as many as 15 generations during the summer ^[25]. In North America, aphids arrive on soybean fields in June, where populations increase by four orders of magnitude, and at the end of the growing season (mid-September), aphids begin the migration back to their overwintering host, reproduce sexually, and overwinter in the egg stage ^[26]. Within Iowa, 40$\%$ of growing seasons from 2004 to 2019 have had populations of aphids large enough to reduce soybean yield, with aphid populations peaking in the middle to the end of August ^[27]. Field surveys in North America have demonstrated that soybean aphid biotypes can co-occur in the same fields ^[28,29]. Laboratory studies have demonstrated that virulent and avirulent biotypes can co-exist on a shared plant for at least 2–3 generations ^[30]. Sub-populations of a herbivorous species can be organized into biotypes, defined as genotypes capable of surviving and reproducing on host plants containing traits (e.g., antibiosis and or antixenosis) conferring resistance to that herbivore ^[31]. Specifically, for the soybean aphid, biotypes are classified based on their ability to colonize soybean varieties expressing Rag-genes (Rag). The (rare) virulent biotype can survive on a soybean host plant containing Rag genes ^[32], while the avirulent (abundant) biotype cannot. On a soybean host plant without Rag genes, the survivability of the avirulent aphid is often greater than the virulent biotype, indicating a clear fitness differential or cost for virulence ^[30].

Global climate change is expected to increase the frequency and intensity of heavy precipitation events in the North-Central United States ^[33], which raises the risk of flooding in soybean crop fields. On these lines, research shows that virulent soybean aphid biotypes have enhanced resilience to other types of stressors, including flooding ^[34] and some pyrethroid insecticides ^[35]. Prolonged flooding events directly damage soybean growth and yield ^[36] and may also alter the quality of soybeans as a host plant for insect herbivores, including changing the plant nutritional quality and/or altering the expression of anti-herbivore defenses. Indeed, soybean aphids experience significant population declines when reared on flood-stress soybeans ^[37]. Host plant flooding stress generally has a negative effect on insect herbivore fitness ^{[38,39,40,41,42]}, although many factors may temper its impact. In soybean aphids, it is possible that resilience to flooding varies between virulent and avirulent soybean aphid biotypes. The mechanisms facilitating virulent aphid colonization of Rag1+2 soybean ^[43] may also confer enhanced resilience to other types of stressors, including flooding; this may introduce a differential fitness response to flooding between aphid biotypes that could alter their population dynamics.

As such, climate change could be a driver of increased differential fitness between avirulent and virulent aphid biotypes on flooded soybean host plants. Scenarios might arise where, on a resistant soybean host plant, the virulent aphid would now be differentially (at least relatively) fitter than its avirulent biotype - that is, the cost of virulence could be much less pronounced. One of our goals is to provide proof of concept models for such scenarios.

2.   Experimental studies and model formulation

Several empirical studies were carried out to test the effect of abiotic stressors on soybean aphids ^[34]. We describe these next.

2.1.   Experimental methods

Laboratory mesocosm assays were conducted to compare how flooding stress in soybeans differentially impacts virulent and avirulent soybean aphid population growth rates. Experiments were conducted using two soybean varieties that were obtained from the National Soybean Research Center (University of Illinois), an aphid-susceptible soybean that contained no Rag genes (LD14-8007) and an aphid-resistant soybean that had pyramided Rag1+Rag2 genes (LD14-8003). Soybean plants were grown from seed in a 24 ℃ environmental growth chamber on a 16.5-hour light: 7.5-hour dark cycle and under optimal water conditions.

At the start of the experiment, 14-day-old plants (V1 growth stage) were placed into individual 1.2 liters (L) of clear plastic jars and randomly assigned to one of two water stress treatments. "Flooded" soybean plants were submerged in $\sim$ 1 L of water, flooding the plant to $\sim$ 1 cm above the soil line, while "Control" plants continued to receive optimal water levels for the duration of the experiment. Plants were then infested with either eight avirulent (biotype 1) or eight virulent (biotype 4) soybean aphid adults and covered with clear, 1 L acrylic tubes to prevent aphid movement between plants.

Aphid density was measured after 10 days by counting the total number of adults and nymphs per plant. Data were subset data by aphid biotype and analyzed with a mixed-model ANOVA using the lme4 and car packages ^[44,45] in R Studio v 4.0.3 (R Core Team). One data point (an avirulent-control-susceptible sample) was excluded from the analysis due to an absence of aphid counts from the sample storage bag; based on similar samples, this error could have been omitted anywhere from 25–100 aphids from the final count. Excluding this observation did not change the analysis results. Each model included soybean variety, water stress, and the soybean-by-water stress interaction as major effects, and the experimental block was included as a random variable. Data were log-transformed to satisfy the assumptions of normalcy and homogeneity of the residuals.

2.2.   Results

Aphid responses to water stress varied between biotypes. Avirulent soybean aphid population densities change significantly in response to water stress treatment (F$_{1, 33.0}$ = 16.38, P $<$ 0.001) and soybean variety (F$_{1, 36.2}$ = 417.5, P $<$ 0.001), although the interaction between water stress and soybean variety was non-significant (F$_{1, 33.0}$ = 0.292, P = 0.593). Indeed, avirulent soybean aphid populations decreased by 28.4 $\%$ with flooding on susceptible soybean and by 44.1$\%$ with flooding on Rag1+2 soybean. Contrasting these results, virulent aphids exhibited no significant response to flooding on either soybean variety (F$_{1, 32.2}$ = 0.275, P = 0.604; Figure 1). Similar to the results in ^[46], virulent aphids also exhibited a slight fitness cost on susceptible soybeans under optimal water control conditions. After 10 days of infestation on control, susceptible soybean, there were 813 ± 124 (mean ± standard error) avirulent aphids per plant compared with 617 ± 88 virulent aphids (a 1.32 ratio of avirulent: virulent). In contrast, this fitness cost was not apparent under flooding conditions, where the ratio of virulent to avirulent aphids was 1.03 (Figure 1).

Based on these experiments, we formulated a competition model for competition among the two aphid biotypes. To this end, we incorporated the differential fitness seen in the experimental studies as an Allee mechanism in one aphid biotype with finite time extinction. We discuss this next.

2.3.   Model formulation

Consider the classical LVM,

where $u$ and $v$ are the population densities of two competing species, $a_1$ and $a_2$ are the intrinsic (per capita) growth rates, $b_1$ and $b_2$ are the intraspecific competition rates, $c_1$ and $c_2$ are the interspecific competition rates. All parameters considered are positive. Note the equilibrium states are achieved only asymptotically, as is the case in many differential equation population models. The dynamics of this system are well studied ^[47]. We recap these briefly,

● $E_0 = (0, 0)$ is always unstable.

● $E_u = (\frac{a_1}{b_1}, 0)$ is globally asymptotically stable if $\dfrac{a_{1}}{a_{2}} > \max\left\lbrace\dfrac{b_{1}}{c_{2}}, \dfrac{c_{1}}{b_{2}}\right\rbrace$. Herein, $u$ is said to competitively exclude $v$ (i.e. $u$ is the strong competitor).

● $E_v = (0, \frac{a_2}{b_2})$ is globally asymptotically stable if $\dfrac{a_{1}}{a_{2}} < \min\left\lbrace\dfrac{b_{1}}{c_{2}}, \dfrac{c_{1}}{b_{2}}\right\rbrace$. Herein, $v$ is said to competitively exclude $u$ (i.e., $v$ is the strong competitor).

● $E^* = \left(\frac{a_1b_2-a_2c_1}{b_1b_2-c_1c_2}, \frac{a_2b_1-a_1c_2}{b_1b_2-c_1c_2}\right)$ exists when $b_1b_2-c_1c_2 \neq 0$. The positivity of the equilibrium holds if $\frac{c_2}{b_1} < \frac{a_2}{a_1} < \frac{b_2}{c_1}$ and is globally asymptotically stable if $b_1b_2-c_1c_2 > 0$. This is said to be the case of weak competition.

● If $b_1b_2-c_1c_2 < 0$, then $E^* = \left(\frac{a_1b_2-a_2c_1}{b_1b_2-c_1c_2}, \frac{a_2b_1-a_1c_2}{b_1b_2-c_1c_2}\right)$ is unstable as a saddle. In this setting, one has an initial condition dependent attraction to either $E_u = \left(\frac{a_1}{b_1}, 0\right)$ or $E_v = \left(0, \frac{a_2}{b_2}\right)$. This is the case of strong competition.

We now provide an approach that incorporates the differential fitness of the two competing biotypes via an Allee effect mechanism. Allee effects are essentially defined as a decline in individual fitness at low population levels that results in extinction once populations fall below a critical level. Several factors, individual and demographic, can cause them, and they can also be driven by predation or human activity. We define certain critical features next ^[48],

Definition 2.1. (Allee threshold) Critical population size or density below which the per capita population growth rate becomes negative.

Definition 2.2. (Component Allee effect) A positive relationship between any measurable component of individual fitness and population size or density.

Definition 2.3. (Demographic Allee effect) A positive relationship between total individual fitness, usually quantified by the per capita population growth rate, and population size or density.

Definition 2.4. (Anthropogenic Allee effect) Mechanism relying on human activity, by which exploitation rates increase with decreasing population size or density: Values associated with the rarity of the exploited species exceed the costs of exploitation at small population sizes or low densities.

Definition 2.5. (Predation-driven Allee effect) A general term for any component Allee effect in survival caused by one or multiple predators whereby the per capita predation-driven mortality rate of prey increases as prey numbers or density decline.

Typical Allee mechanisms, in the single species case, would be of the form,

where $A$ is the Allee threshold; if $A > 0$, there is a strong Allee effect, and if $A \leq 0$, there is a weak Allee effect - this means that the growth rate is low at low population levels, but there is no Allee threshold. The Allee effect has large-scale consequences in population management ^[48] and has been considered in mating systems ^[49] and in biocontrol applications ^[49,50]. The Allee effect in competitive systems began with the seminal work of Wang et al. ^[51] who considered the Lotka-Volterra population model, with the Allee effect modeled as,

where $C = 0$ is no Allee effect, and the greater the $C$, the greater the effect. In this setting, the Allee effect is in weak form. Various rich dynamics were discovered, including the existence of multiple interior equilibrium points. They later extended their results to the metapopulation systems ^[52]. The Allee effects' interplay with dispersal in competition systems has also been investigated ^[53].

We have also used the idea of a strong Allee effect, with the Allee type threshold, to model avirulent and virulent aphids on a resistant/Rag-type soybean plant ^[54]. However, there are several problematic issues with this approach.

● The threshold value of $A = 30$, used in ^[54], and determined in earlier experiments ^[55] cannot predict the current experimental data. Earlier experimental results predict that avirulent aphids arriving in initial numbers under 30 would not be able to colonize a resistant variety via feeding facilitation alone ^[55]. However, this is not seen in current experiments; see Figure 11(b).

● Thus, the exact value of the threshold $A$ is unknown, ^[32].

● Based on current and past empirical studies ^[32,34], we require a strong type effect, with density dependence ^[34]–but setting $A = A_{0}$ (where $A_{0}$ is decided ad-hoc), does not provide this.

● Below this threshold, the populations should be driven to extinction rapidly (such as seen with avirulent populations in low numbers on a rag-type plant ^[32]). However, the strong effect used earlier ^[54] enables only for an exponential decay in the populations, while the population decrease actually occurs very quickly, so it is more in tune with finite time extinction ^[56].

To the best of our knowledge, there is no consideration of an Allee effect in competitive systems with a finite time extinction mechanism (FTEM). To this end, we define an Allee mechanism in a single species setting via,

Here, $0 < q \leq 1, 0 \leq p < 1$. The parameter $q$ describes the intensity of an environmental effect (for example, flooding) due to a changing climate. $q = 1$, is the no flood case, with $q \approx 0$ being the strongest possible flood case. Here, there would be no growth rate essentially, and populations would decay $\approx -c v^{p}$. The reason being in this setting ($q \approx 0$), (2.4) reduces to, $\frac{dv}{dt} = - b_{2}v^{2} - c v^{p} \leq - c v^{p}$. This follows via a simple comparison argument. Thus, the term $- c v^{p}$ is the dominant term in the reaction kinetics, and solving $\frac{dv}{dt} = -c v^{p}$ yields finite time extinction of the $v$ population. This can be proven very similar to Lemma 3.5, and thus the $v$ population decays $\approx -c v^{p}$.

One of our chief motivations is that flooding, apart from having an Allee effect with FTEM, also reduces the growth rate of the population; hence, the term $q$ multiplies the growth rate. The parameters $c$ and $p$ control the density-dependent Allee threshold regarding plant resistance and biotype. In theory, increasing $c$ causes more initial data to be attracted to extinction equilibrium, and increasing $p$ increases the rate of attraction. Thus, $c, p$ would be higher on a resistant plant for the avirulent biotype than a virulent biotype. However, on a susceptible plant $c, p$ would be higher for a virulent aphid biotype than an avirulent one due to the cost for virulence ^[32]. For tractability and to reduce parameters, one could choose $c = 1$ and just vary $p$ as well.

Thus, the Allee threshold model with FTEM we propose is driven by

● Current empirical studies ^[34].

● The effect of a changing climate ^[33].

● Plant virulence ^[55].

Motivated by the above, we now present a Lotka Volterra competition model. The fitness dynamics of $u$ (virulent biotype) are governed as earlier. The Allee mechanism models a fitness differential between the competing biotypes, where the Allee mechanism is present only in the $v$ population. This is seen in empirical evidence ^[34] as the avirulent biotype is strongly affected by flooding– while no significant difference is found in the virulent biotype population, on flood vs. control setting (see Figure 1). Our model is as follows,

where $0 < p < 1$, $0\leq q \leq 1$.

Remark 1. Since $q = 1$ provides a classical Lotka-Volterra competition model, we focus our analysis on $q\in(0, 1)$.

2.4.   Synopsis of current manuscript

The following is accomplished in this manuscript,

● A two-species ODE Lotka-Volterra competition model, with a non-smooth degradation term, to model differential fitness via an Allee effect in two competing sub-populations is introduced. This model exhibits finite time extinction mechanisms (FTEM).

● This model can lead to the dynamics of bi-stability where the only interior equilibrium point is a saddle, which can be seen in Figure 6 and later proved in Lemma 3.9 and Theorem 3.10.

● FTEM can also enable a competitor to avoid competitive exclusion and persist for various regimes of initial conditions. This is later seen via Theorem 3.4 and can be observed in Figure 2(a). To this end, the differential fitness needs to affect only a portion of the weaker competitors' population, which is proven in Theorem 3.4.

● Rich dynamical and bifurcation phenomena are uncovered as can be seen in Figures 7 and 8 and worked on in Section 3.

● The requirement on the smallness of initial data to yield FTE is described via Theorem 3.6. Also, see Figure 5.

● Recent observations and experimental evidence are described and discussed, showing that there is a fitness differential between virulent and avirulent biotypes driven by abiotic stressors such as flooding events, which can be graphically observed in Figure 1.

● A mathematical framework is developed to model scenarios of competing aphid biotypes in Section 6, based on the experimental evidence. A separate model is proposed and simulated, where novel dynamics are observed and portrayed in Figures 11 and 12

● We discuss our results in the context of managing the invasive pestiferous soybean aphid under climate change.

3.   Dynamical analysis

3.1.   Existence of equilibrium points

In this subsection, we analyze the existence of equilibrium of system (2.5). We are interested in studying the dynamics of system (2.5) in the closed first quadrant $\mathbb{R}^2_+$ in the (u, v) plane, keeping biological implications for the system in mind. It is obvious that $E_0(0, 0)$ is an equilibrium point of system (2.5). We can have two types of boundary equilibrium points: $E_u\left(\frac{a_1}{b_1}, 0\right)$, which is trivial, and $E_v(0, \bar{v})$ as described in Lemma 3.1. The whole plane $\mathbb{R}_+^2$ is the positive invariant set under parameter $q$ for system (2.5).

In this section, we discuss the existence of interior equilibrium points of system (2.5). In order to get the equilibrium points of system (2.5), we consider the nullclines of both the population $u$ and $v$ of system (2.5), which are given by:

Simplifying system of equations (3.1) to get an explicit form in terms of population $v$, we get,

We study the dynamics of polynomial (3.2) to find the number of equilibrium points possible for system (2.5).

Lemma 3.1. Let $v_{\phi}^{(p-2)} = \dfrac{b_2}{c(1-q)(1-p)}$. Then we have:

1) If $q < \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there does not exist any boundary equilibrium of the form $E_v(0, \bar{v})$.

2) If $q = \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists an unique boundary equilibrium $E_v(0, \bar{v})$.

3) If $q > \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists two boundary equilibrium $E_{v_1}(0, \bar{v_1})$ and $E_{v_2}(0, \bar{v_2})$.

Proof. The detailed proof can be found in Appendix (A.1) □

Lemma 3.2. If $b_1b_2-c_1c_2 \leq 0$ (strong competition), there exists only one unique interior coexistence equilibrium.

Proof. The detailed proof can be found in Appendix (A.2) □

Lemma 3.3. When $b_1b_2-c_1c_2 > 0$ (weak competition) and $v_{max}^{p-2} = \dfrac{b_1b_2-c_1c_2}{c(1-q)(1-p)}$ then we have:

1) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} < 0$, then there does not exist any positive interior equilibrium

2) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} = 0$, then there exists one unique interior equilibrium

3) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} > 0$, then there exists two positive interior equilibrium points.

Proof. The simplification of the nullclines gives

We attempt to study the dynamics of the slope of $\phi(v)$ to determine the existence of equilibrium points. We have,

As $(b_1b_2-c_1c_2) > 0$, $\phi'(v)$ can be of any sign depending on the magnitude of $v$. Let us compute an extremum $v_{max}$ using $\phi'(v_{max}) = 0$. Solving for the extremum we get, $v_{max}^{p-2} = \frac{b_1b_2-c_1c_2}{cb_1(1-q)(1-p)}$. As $b_1b_2-c_1c_2 > 0$ and $0 < p, q < 1$ thus $v_{max}$ is positive.

To study the nature of the extrema, we find the second derivative, which gives:

As $0 < p, q < 1$ so $\phi''(v_{max})$ is strictly negative for any value of $p$ and $q$. Thus, at $v_{max}$, we have a maximum. Now, as the parameters are fixed so, $v_{max}$ is a unique positive root. So, $v_{max}$ is a global maximum. Then, the maximum value of the function is

We see $\phi(v) \to -\infty$ as $v\to 0$ and $v\to \infty$. So the curve is downward facing as seen in Figure 2.

Thus, if $\phi(v_{max}) < 0$, then it is obvious that there does not exist any positive root of $\phi(v)$. If $\phi(v_{max}) = 0$ then $v_{max}$ is the only root of $\phi(v)$. As $v_{max}$ is the only maxima if $\phi(v_{max}) > 0$, then we can have two distinct positive roots of $\phi(v)$. Hence, the number of equilibrium points is determined for the polynomial $\phi(v)$ when $b_1b_2-c_1c_2 > 0$. □

Via Lemmas 3.1–3.3, we have the existence of equilibrium points for system (2.5). This can be summarized via Theorem 3.4 as:

Theorem 3.4.

1) System (2.5) has a trivial equilibrium $E_0(0, 0)$.

2) System (2.5) has a boundary equilibrium $E_u(a_1/b_1, 0)$ where species $u$ persists.

3) System (2.5) has $E_v(0, \bar{v})$ as boundary equilibrium where species $v$ persists.

When $v_{\phi}^{p-2} = \dfrac{b_2}{c(1-q)(1-p)}$, the number of equilibrium of the form $E_v(0, \bar{v})$ can be determined by:

(a) If $q < \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there does not exist any boundary equilibrium of the form $E_v(0, \bar{v})$.

(b) If $q = \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists an unique boundary equilibrium $E_v(0, \bar{v})$.

(c) If $q > \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists two boundary equilibrium $E_{v_1}(0, \bar{v_1})$ and $E_{v_2}(0, \bar{v_2})$.

4) If $(b_1b_2-c_1c_2) \leq 0$, then system (2.5) has an unique interior equilibrium.

5) If $(b_1b_2-c_1c_2) > 0$ and $v_{max}^{p-2} = \dfrac{b_1b_2-c_1c_2}{c(1-q)(1-p)}$, then

(a) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} < 0$, then system (2.5) has no equilibrium as in Figure 2(a).

(b) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} = 0$, then system (2.5) has an unique interior equilibrium as in Figure 2(b).

(c) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} > 0$, then system (2.5) has two different interior equilibrium points as in Figure 2(c).

Remark 2. System (2.5) shows due to competition, there can be extinction of one of the species while the other species persists. The population without the effect of flooding $v$ can go extinct in finite time, as can be seen in Figure 3 and 4. The coexistence of both species is also possible due to the existence of interior equilibrium points, as can be seen in Figure 2.

3.2.   Stability analysis of the equilibrium points

In this subsection, we study the dynamics of system (2.5) in the neighborhood of each equilibrium. We study the stability of each equilibrium using the Jacobian matrix $J(u, v)$ of system (2.5) given by,

$J(u, v) = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} = \begin{bmatrix} a_1-2b_1u-c_1v & -c_1u \\ -c_2v & a_2q-2b_2v-c_2u-cp(1-q)v^{(p-1)} \end{bmatrix}$

As the system has a singularity at $v = 0$, certain stability analysis of boundary equilibrium points cannot be performed using classical techniques. Thus, the stability analysis herein is mostly shown using calculus techniques by graphical interpretation.

$E_0(0, 0)$ is unstable. However, showing this via standard linear analysis is not possible due to the singularity in the $v$ component in $a_{22}$ in the Jacobian matrix above. However, this is seen, as plugging in $u = 0$ in (2.5) will yield finite time extinction in the $v$ component. However, we know there is data that is attracted to $(\frac{a_{1}}{b_{1}}, 0)$, in finite time via Theorem 3.6, and then $u \rightarrow \frac{a_{1}}{b_{1}}$ as $t \rightarrow \infty$. Thus, the $v$-axis is the stable manifold for $(0, 0)$, while the $u$-axis is the unstable manifold–thus, $E_0(0, 0)$ is unstable as a saddle. However, the stable manifold is not unique due to the lack of smoothness of the system ^[57].

Next, we investigate the stability of the boundary and the interior equilibrium points. The analysis of the boundary equilibrium $E_u(\frac{a_1}{b_1}, 0)$ is not possible via linear stability analysis due to the lack of differentiability of the system (2.5) at $v = 0$. Nonetheless, the following result can be provided:

Lemma 3.5. Consider the boundary equilibrium $E_u\left(\frac{a_1}{b_1}, 0\right)$. This is locally stable. That is, there exists certain data that is attracted to this equilibrium in finite time.

Proof. Consider the equation for $v$,

via positivity and simple comparison, we have,

We consider the auxiliary equation,

Note via simple comparison, $v \leq V$, for all time $t$. Set $V = \frac{1}{U}$, then,

This yields,

Clearly, via comparison with $\frac{dW}{dt} = c(1-q)W^{2-p} - a_2 q W, \ 0 < p < 1$, we see that $U$ will blow up in finite time for sufficiently large initial condition $U(0)$. That is if $U(0)$ is chosen, such that

Then

However, since $V = \frac{1}{U}$, we have,

That is, $V$ will go extinct in finite time, for initial data small enough, that is given by,

Since $v\leq V$, $v$ will go extinct in finite time as well. That is,

For initial condition chosen such that

Now, the $u$ equation reduces to,

which implies,

This proves the lemma. □

We next give an explicit characterization of the initial conditions that yield finite time extinction.

Theorem 3.6. Consider positive initial conditions $(u_{0}, v_{0})$ s.t., $u_{0} < \frac{a_{2} q}{c_{2}} < \frac{a_{1}}{b_{1}}$, $v_{0} < \frac{a_{2} q}{b_{2}}$, $\left(\frac{a_{2}q}{b_{2}}\right)^{(1-p)} c_{2} u_{0} \left(e^{\frac{c_{1} v_{0}}{a_{2}q}} \right) < c(1-q)$ and $(a_{2} q) (v_{0})^{1-p} < \left(c(1-q)(1-p) - (1-p)\left(\frac{a_{2}q}{b_{2}}\right)^{(1-p)} c_{2} u_{0} \left(e^{\frac{c_{1} v_{0}}{a_{2}q}} \right)\right)$, then for any solutions initiating from these initial conditions, $v$ goes extinct in finite time, that is,

Proof. Consider the equation for $u$,

If initial conditions are chosen s.t. $u_{0} < \frac{a_{1}}{b_{1}}$, then by simple comparison with the logistic equation, $u < \frac{a_{1}}{b_{1}}$, for all time $t$. Thus, $a_{1}u - b_1 u^2 \geq 0$, and we have,

This follows as by simple comparison, $v \leq v_{0}e^{a_{2}qt}$, for all $t$. Thus we have,

Now, via the positivity of solutions from the $v$ equation, we can ascertain,

dividing by $v^p$, yields,

Thus, we obtain,

Now, if initial data and parameters are chosen such that $u_{0} < \frac{a_{2} q}{c_{2}} < \frac{a_{1}}{b_{1}}$, then we have $\left(a_{2}q - c_{2} u \right) > 0$, for all time $t$. Thus, using the lower bound on $u$, we have,

via comparison, and the upper bound on $v$ via logistic control, that is $v \leq \frac{a_{2}q}{b_{2}}$, for all time as long as $v_{0} \leq \frac{a_{2}q}{b_{2}}$, it follows that,

Thus,

Integration of the above inequality in time yields,

Now, standard theory yields finite time extinction of $v$ if initial conditions are chosen such that, $c(1-q) > \left(\frac{a_{2}q}{b_{2}}\right)^{(1-p)} c_{2} u_{0} \left(e^{\frac{c_{1} v_{0}}{a_{2}q}} \right)$ and,

Then $v$ will go extinct in finite time. This proves the theorem. □

Remark 3. Theorem 3.6 aims to show that there are certain initial populations, from which if we initiate our competing population system (2.5), then the $v$ population in particular will go extinct in a finite time. To enable this, certain parametric restrictions need to be enforced, restricting the initial populations in terms of the model parameters. In essence, we need the $u$ population to be strictly bounded above for all time and the $v$ population to start fairly small. The mechanics of the theorem require us to enforce $a_{1}u - b_1 u^2 \geq 0$, or $u \leq \frac{a_{1}}{b_{1}}$. In order for the $u$ population to be strictly bounded above by this carrying capacity for all time, we need to start it lower than its carrying capacity, which is $\frac{a_{1}}{b_{1}}$. In this case, via comparison with the logistic equation, we can enforce that the $u$ population will remain below $\frac{a_{1}}{b_{1}}$ for all time. Without this assumption, even if $u$ decays to a steady state below the critical $\frac{a_{1}}{b_{1}}$ level, one is not able to ascertain an explicit bound, below which we can maintain the $u$ population for all time. Thus, this is an essential hypothesis in the theorem.

Lemma 3.7. Let $v_{\phi}^{p-2} = \dfrac{b_2}{c(1-q)(1-p)}$ then

1) If $q = \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists an unique boundary equilibrium $E_v(0, v_\phi)$ and it is a saddle.

2) If $q > \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists two boundary equilibrium $E_{v_1}(0, v_1)$ and $E_{v_2}(0, v_2)$ and if $\frac{a_1}{c_1} < v_1 < v_2$, then $E_{v_1}$ is a saddle while $E_{v_2}$ is a stable node.

Proof. The proof follows in Appendix (A.3) □

Lemma 3.8. If $(b_1b_2-c_1c_2)\leq 0$, then the unique interior point is always an unstable equilibrium point.

Proof. The detailed proof can be found in Appendix (A.4) □

Lemma 3.9. If $b_1b_2-c_1c_2 > 0$ and $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} > 0$ then there exists two positive interior point ($E_1 < E^*$) where $E^*(u^*, v^*)$ is a stable if q is in the range,

where $v_{max}^{p-2} = \dfrac{b_1b_2-c_1c_2}{c(1-q)(1-p)}$ else it is a saddle point.

Proof. The proof follows in Appendix (A.5) □

According to Lemmas 3.7, 3.8, and 3.9, we have an understanding of the stability of the equilibrium points of the system (2.5). The results can be summarized into a theorem given by:

Theorem 3.10. The study of the stability of the system (2.5) shows,

1) The trivial equilibrium $E_0(0, 0)$ is always unstable.

2) The boundary equilibrium $E_u\left(\frac{a_1}{b_1}, 0\right)$ is stable when there exists an interior saddle $E(u^{*}, v^{*})$. Initial data below $W^{s}(E)$, the stable manifold of this saddle, is attracted to $E_u\left(\frac{a_1}{b_1}, 0\right)$ in finite time.

3) The boundary equilibrium $E_v(0, \bar{v})$ is stable when $v_{\phi}^{p-2} = \dfrac{b_2}{c(1-q)(1-p)}$, so

(a) If $q = \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists an unique boundary equilibrium $E_v(0, v_\phi)$ and it is a saddle point.

(b) If $q > \dfrac{b_2v_\phi(2-p)}{a_2(1-p)}$, then there exists two boundary equilibrium $E_{v_1}(0, v_1)$ and $E_{v_2}(0, v_2)$ and if $\frac{a_1}{c_1} < v_1 < v_2$, then $E_{v_1}$ is a saddle while $E_{v_2}$ is a stable node.

4) When $b_1b_2-c_1c_2\leq 0$, then the unique interior point is always an unstable point.

5) When $b_1b_2-c_1c_2 > 0$ and $v_{max}^{p-2} = \dfrac{b_1b_2-c_1c_2}{c(1-q)(1-p)}$, then there exists two interior equilibrium point ($E^*_1(u_1, v_1) < E^*_2(u^*, v^*)$),

(a) If $q < \dfrac{1}{a_2b_1(1-p)}\left\{(b_1b_2-c_1c_2)(2-p)v_{max}+c_2a_1(1-p)\right\}$, then interior equilibrium points dont exist.

(b) If $q = \dfrac{1}{a_2b_1(1-p)}\left\{(b_1b_2-c_1c_2)(2-p)v_{max}+c_2a_1(1-p)\right\}$, then the equilibrium points undergoes Saddle Node Bifurcation.

(c) If $(a_2b_1q-c_2a_1)-v_{max}(b_1b_2-c_1c_2)\dfrac{(2-p)}{(1-p)} > 0$, then the equilibrium point $E^*_1$ is a saddle and $E_2^*$ is stable if following conditions hold

Remark 4. System (2.5) shows that complete extinction is not possible while the persistence of species $u$ and extinction of $v$ is possible based on the initial population threshold. The competitive exclusion of population $u$ by species $v$ is possible and depends on the amount of flooding. We can have almost two interior equilibrium points. According to Lemma 3.9, we see that when $E_{max}(u_{max}, v_{max})$ is the equilibrium and we cross from one side of a plane to another by change of parameter $q$ the system (2.5) changes from no interior equilibrium to two equilibrium points, which can be seen in Figure 7 and Figure 8(a). Thus, there can exist a saddle-node bifurcation at $E_{max}(u_{max}, v_{max})$, which we study in the following section.

3.3.   Saddle-node bifurcation

When $(b_1b_2-c_1c_2) > 0$, then by Lemma 3.3, there can be two interior equilibrium points from no equilibrium point when the nullclines crosses $E_{max}(u_{max}, v_{max})$ equilibrium point due to change in parameter $q$.

The Jacobian matrix for such an equilibrium point $E_{max}(u_{max}, v_{max})$ is,

Thus,

As $v_{max}$ and $u_{max}$ are not zero then,

It is obvious that $\lambda_1 = 0$ and $\lambda_2 = -b_1u_{max}-v_{max}(b_2-(b_1b_2-c_1c_2)/b_1)$ are the eigenvalues of the $J(E_{max})$. Thus, if $-a_1+v_{max}(c_1-b_2+(b_1b_2-c_1c_2)/b_1) \neq 0$, then $\lambda_2 \neq 0.$ Thus, from the conditions, we can obtain that,

is a saddle-node bifurcation surface.

Sotomayor's theorem ^[58] is used to verify the transversality conditions for the occurrence of saddle-node bifurcation of the parameter on the surface $SN_1$. We know that $J(E_{max})$ has one simple zero eigenvalue. If V and W represent the eigenvectors for the zero eigenvalues of the matrix $J(E_{max})$ and $J(E_{max})^T$, respectively, then V and W are,

Obviously, V and W satisfy the transversality conditions,

$W^TF_q(E_{max}; SN_1) = \dfrac{-b_1u_{max}}{c_1}(a_2+p) \neq 0$,

$W^TD^2F(E_{max}; SN_1)(V, V) = \dfrac{-2b_1^2u_{max}}{c_1^2c_2v_{max}}(b_1b_2-c_1c_2)(p-2) \neq 0$, as $b_1b_2-c_1c_2 > 0$ and $0 < p < 1$.

For the occurrence of the saddle-node bifurcation at the parameters on the surface $SN_1$. Hence, it can be stated that when the parameters cross from one side of the surface to the other side, the number of equilibrium points in the model changes from zero to two, and the two equilibrium points, which are interior equilibrium points, are hyperbolic saddle and node.

The saddle-node bifurcation existence can be numerically seen for a particular parameter set as in Figure 7.

3.4.   Pitchfork bifurcation

Theorem 3.11. Consider system (2.5), when the boundary equilibrium points $E_v(0, \bar{v})$ satisfies the nullcline equation $qa_2-b_2\bar{v}-c(1-q)\bar{v}^{(p-1)} = 0$ then a pitchfork bifurcation occurs as $q \to q^*$, where,

Proof. We obtain the bifurcation parameter by studying the gradient of the nullclines.

The nullclines from system (2.5) are,

The slope of the nullclines is to be determined at the equilibrium points $E_v(0, \bar{v})$. The explicit form of $\bar{v}$ can be derived from the v-nullclines as it is a root of the equation $qa_2-b_2\bar{v}-c(1-q)\bar{v}^{(p-1)} = 0$. The gradient of the nullclines at $E_v(0, \bar{v})$ is given by $\dfrac{df}{dv}\rvert_{v = \bar{v}} = -\dfrac{c_1}{b_1}$ and $\dfrac{dg}{dv}\rvert_{v = \bar{v}} = -\dfrac{b_2}{c_2}+\dfrac{c(1-q)(1-p)}{c_2}\bar{v}^{(p-2)}$. When pitchfork bifurcation occurs, the interior equilibrium points collide with the boundary equilibrium, and the slope of the nullclines are equal. Thus, $-\dfrac{c_1}{b_1} = -\dfrac{b_2}{c_2}+\dfrac{c(1-q)(1-p)}{c_2}\bar{v}^{p-2}$. Simplifying the equality with the fact that $qa_2-b_2\bar{v}-c(1-q)\bar{v}^{(p-1)} = 0$ yields the following result,

As $q$ decreases to $q^*$, the $g(v)$ nullcline moves downward, and the two interior equilibrium points and one boundary equilibrium come closer together. The shape of the nullclines follows this under the parametric restriction enforced. Since the two slopes of prey and predator nullclines at x = 0 are chosen to be the same, by continuity, the three equilibrium points merge into one as $q \to q^*$. Now, as $q$ is further decreased, the $g(v)$ nullcline is completely below the prey nullcline, $g(x) < f(x), \forall x$, and there is a single boundary equilibrium. Thus, by definition, pitchfork bifurcation has occurred. For a set parameter space, pitchfork bifurcation can be seen in Figure 8(b). □

4.   Comparison with the classical competition model

System (2.5) can be trivially transformed to the LVM 2.1 if $q = 1$. We study different cases of how the dynamics of (2.5) differ from the LVM by changing parameters $p$ and $q$ in the range $(0, 1)$.

4.1.   Competitive exclusion

LVM yields that when $\dfrac{a_{1}}{a_{2}} < \min\left\lbrace\dfrac{b_{1}}{c_{2}}, \dfrac{c_{1}}{b_{2}}\right\rbrace$, then $E_v = \left(0, \dfrac{a_2}{b_2}\right)$ is globally asymptotically stable. Thus, for any parameter set in the above regime with $q = 1$, $v$ competitively excludes $u$.

If we choose a $q$ very close to 1, then the dynamics can change. The parameters chosen to satisfy the competitive exclusion condition can be seen in Figure 9(a), which yields $E_v = (0, \frac{a_2}{b_2})$ to be globally asymptotically stable (with $q = 1$). With the choice of $q = 0.91$ and $p = 0.6$, we see that $E_v = (0, \frac{a_2}{b_2})$ becomes unstable, and we get two stable equilibrium points, i.e., $(u^{*}, 0)$ boundary equilibrium and an interior equilibrium, on the two sides of the stable manifold of the interior saddle, as seen in Figure 9(b). Thus, for any initial data below the stable manifold, $(u^{*}, 0)$ boundary equilibrium is stable whereas the $v$ population goes to extinction, while above the stable manifold, coexistence takes place where the $E^*$ equilibrium becomes stable. Thus, the $u$ population is always persistent in the system and does not go extinct, which is completely different from the case in the classical LVM.

4.2.   Weak competition

The classical LVM also yields weak competition type dynamics under parametric restrictions $\frac{c_2}{b_1} < \frac{a_2}{a_1} < \frac{b_2}{c_1}$ and $b_1b_2-c_1c_2 > 0$, wherein the coexistence of two species is possible and globally asymptotically stable over time. In the LVM model (2.1), we can get one interior equilibrium that can be asymptotically stable under these restrictions, which in system 2.5 can be achieved by satisfying these parametric restrictions and with $q = 1$ as seen in Figure 9(c).

In system 2.5, if we use $0 < q < 1$, we can get two interior equilibrium points with one locally stable interior equilibrium point and the other as an interior saddle equilibrium point. Due to the saddle dynamics, we can have bi-stability where the stable manifold of the saddle equilibrium divides the phase space into two regions; depending on the initial condition, the populations can either be attracted to the interior equilibrium or can be attracted to the $(u^{*}, 0)$ boundary equilibrium as seen in Figure 9(d).

5.   The spatially explicit/PDE case

5.1.   Preliminaries

The following results are classical, (see ^[59] for details).

Lemma 5.1. Consider a $2\times 2$ system of reaction-diffusion equations, where each equation is defined as follows: For all $i = 1, 2$

where $d_i \in (0, +\infty)$, $f = (f_1, f_2):\mathbb{R}^2 \rightarrow \mathbb{R}^2$ is continuously differentiable on $\Omega$ and $u_{i0}\in L^{\infty}(\Omega)$. Then, there exists a time interval $T > 0$ within which a unique classical solution to (5.1) exists, i.e., the solution is well-defined and smooth on $[0, T)$. Let $T^*$ be the maximum value for all such intervals $T$. It follows that

If the non-linearity $(f_1, f_2)$ is also quasi-positive, i.e., if

then

Lemma 5.2. Under the same notation and assumptions as in Lemma 5.1, let us consider an additional condition. Suppose that $f$ exhibits at most polynomial growth and there exist $\mathbf{b}\in \mathbb{R}^m$ and a lower triangular invertible matrix $P$ with non-negative entries such that for any $r \in [0, +\infty)^m,$ we have

Then, for any initial value $u_0 \in L^{\infty}(\Omega, \mathbb{R}_+^m),$ the system 5.1 admits a strong global solution.

Based on the given assumptions, it is widely recognized that the following local existence result, due to Henry ^[59], holds true:

Theorem 5.3. The system (5.1) possesses a unique and classical solution $(u, v)$ defined over the interval $[0, T_{\max}]\times \Omega$. If $T_{\max} < \infty$, then

where $T_{\max }$ denotes the eventual blow-up time in $\mathbb{L}^{\infty}(\Omega).$

5.2.   The constant growth coefficient case and Turing instability

In this section, we study the spatially explicit PDE versions of systems (2.5), where the coefficients are constants. Note first, the system under consideration,

where $\nabla u \cdot n = \nabla v \cdot n = 0$, $u(x, 0) = u_{0}(x), v(x, 0) = v_{0}(x)$, $0 < p < 1$, $0\leq q \leq 1$.

Definition 5.4. A solution to (5.3), is a non-negative classical solution $(u, v)$ with $u, v \in C([0, T^{*}) \times \overline{\Omega}) \cap C^{1, 2}((0, T^{*}) \times \overline{\Omega})$, where $T^{*} \in (0, \infty]$ denotes the maximal existence time. Moreover, for any solution $(u, v)$ if $T^{*} < \infty$, then $(u, v)$ is said to blow up in finite time, that is,

Remark 5. It is known that for systems of type (5.3), there exist classical solutions locally, which could be extended to global solutions. This follows via standard application of Lemmas 5.1–5.2. However, in the non-Lipschitz case of $0 < p < 1$, uniqueness may or may not hold depending on parameter restrictions. This follows via the application of the techniques in ^[60].

We next state the following result concerning (5.3),

Lemma 5.5. Consider the reaction-diffusion system (5.3). There is no Turing instability in the system in case of no flooding, i.e., $q = 1$.

Proof. A standard linearization of (5.3) yields,

$J_{22} < 0$ and $J_{11} < 0$ as $u^*, v^* > 0$, hence $J_{11}J_{22} > 0$ and so the necessary condition required for Turing instability (that is $J_{11}J_{22} < 0$), is violated. This proves the Lemma. □

We first state the following result.

Suppose ($u^*, v^*$) is a stable interior equilibrium point. We investigate the possibility of destabilizing this stable interior equilibrium point via unequal diffusion. This will demonstrate the existence of Turing instability caused by flooding-mediated competition, which is clearly not possible without competition via Lemma 5.5.

The system (2.5) can be rewritten as

where $X = [u\; \; v]$ and D = diag($d_1, d_2$) and the Jacobian $J$ is defined in (5.4).

We analyze the system with Neumann conditions on the boundaries on a spatial domain $\Omega = [0, L] \times [0, L]$ and the wave number $k = n\pi/L$. Since ($u^*, v^*$) is a stable interior equilibrium point, when there is no diffusion present i.e., $D = 0$, we have the following condition: $Trace(J) = J_{11}+J_{22} < 0$ and $\det(J) = J_{11} J_{22}-J_{12} J_{21} > 0$.

The solution to the system (5.5) is of the form

where $c_k$ are amplitudes determined by Fourier expansion and $X_k(\Omega)$ are the eigenvectors dependent on the eigenvalues $\lambda(k)$ of the eigenvalue problem where $k = n\pi/L$. Thus, to obtain Turing instability, we need at least one positive eigenvalue of the system (5.3)

where $\lambda$ denotes the eigenvalue. For simplicity, we normalize the diffusion coefficients to obtain D = diag(1, d). From condition (5.7), we obtain that the eigenvalues are the roots of the equation

where

Now, to get a positive eigenvalue which shows instability due to diffusion, the following must hold

from (5.10) and the fact that $Trace(J) < 0$, we can conclude that we need $J_{11}J_{22} < 0$.

and for a fixed value of $d$ the value of $k$ can be determined by $h(k^2) < 0$. Thus, we state the following lemma,

Lemma 5.6. Consider the reaction-diffusion system (5.3), when $c > 0$. In the case of an Allee effect, i.e., $q\in (0, 1)$, the necessary condition for Turing instability is if $b_1 > 1$ and $q > q^*$ where $q^* = \frac{1}{a_2b_1(1-p)}\left(v^*b_2+(b_1b_2-c_1c_2)(1-p)v^*+c_2a_1(1-p)\right)$.

Proof. The Jacobian matrix used is

The necessary condition for Turning instability is $J_{11}J_{22} < 0$. As $J_{11} < 0$, thus, to have Turning instability, we need $J_{22}$ to be positive.

Thus, simplifying $J_{22} > 0$

Using the nullcline equation,

Thus,

In order to have Turing instability, Eq (2.5) stability conditions must hold.

From Theorem 3.10, we have

Thus, $\bar{q} > q > q^*$.

Simplifying, we have $b_1 > 1$. This proves the lemma. □

Remark 6. The condition $b_1 > 1$ is biologically not feasible as the rate of intraspecies competition lies in the region $[0, 1]$. However, the conditions state that diffusion can create instability in the system under some level of flooding.

5.3.   The nonconstant growth coefficient case

The spatially inhomogeneous problem has been intensely investigated in the past 2 decades ^{[1,61,62,63,64]}. The premise here is that $u, v$ do not have resources that are uniformly distributed in space; rather, there is a spatially dependent resource function $m(x)$. We consider again a normalized generalization of the classical formulation,

Note that $q = 1$ is the classical case. We consider $m$ to be non-negative on $\Omega$ and bounded. We recap a seminal classical result ^[61], which shows that the slower diffuser wins,

Theorem 5.7. Consider (5.13) and (5.14), when $q = 1$, and $d_{1} > d_{2}$, solutions initiating from any positive initial data $(u_{0}(x), v_{0}(x))$ converge uniformly to $(0, v^{*}(x))$.

That is, the slower diffuser wins in the case of equal kinetics. However, when $0 < p, q < 1$, that is, the Allee effect is active in the competitive system as a finite-time extinction mechanism, the slower diffuser can lose, depending on the initial conditions. In the ensuing analysis the value of the constant $"C"$ can change from line to line, and sometimes within the same line if so required. Also, the existence of classical solutions to (5.13) and (5.14) follows via standard application of Lemmas 5.1 and 5.2, as well as the techniques in ^[60]. We next state the following result,

Theorem 5.8. Consider (5.13) and (5.14), where $\Omega \subset \mathbb{R}^{n}$, $n = 1, 2$, is a bounded domain, with smooth boundary, $q = 1$, and $d_{1} > d_{2}$. There exists positive initial data $(u_{0}(x), v_{0}(x))$, for which the solutions converge to $(0, v^{*}(x))$. However, any $0 < q < 1$ solutions with the same diffusion coefficients, initiating from the same data, will converge to $(u^{*}(x), 0)$ in finite time for a sufficiently chosen $p \in (0, 1)$.

Proof. Via comparison with the logistic equation ^[1], we see that $u \leq C m(x)$, $\forall x, t \in \Omega \times [0, \infty)$.

We now multiply the $v$ equation in (5.13) by $v$ and integrate by parts to obtain

Using positivity of $u, v$ we obtain

Then, it follows that,

Choosing $C_{1} = \min{(d_{1}, c(1-q))}$ will yield,

Our next goal is to show that

for some $0 < \alpha < 1$, as this would yield,

Then, we will have the finite time extinction of $v$ in comparison to the ODE

Now, recall the Gagliardo-Nirenberg-Sobolev (GNS) inequality ^[65],

for $\phi \in W^{m, q^{'}}(\Omega)$ provided $p^{'}, q^{'}, q > 1, 0 \leq \theta \leq 1$, and

Now, consider exponents s.t.

for $0 < p < 1$. This yields

as long as

and

We raise both sides of (5.25) to the power of $l$, $0 < l < 2$, to obtain

Using Young's inequality on the right-hand side (for $a b \leq \frac{a^{r}}{r} + \frac{b^{m}}{m}$), with $r = \frac{2}{l \theta}, \ m = \frac{q^{*}}{l (1-\theta)}$, yields

We notice that given any $1 < q^{*} < 2$, it is always possible to choose $0 < l < 2$, and $0 \leq \theta \leq 1$ s.t, $\frac{1}{r} + \frac{1}{m} = 1$,

by choosing

thus for any given $1 < q^{*} < 2$, we need to choose $l$ s.t,

or

while being mindful that $0 \leq \theta \leq 1$. This is for $n = 1$. For $n = 2$ case, given $1 < q^{*} < 2$, we need to choose $l$ s.t,

or

Whilst again being mindful that $0 \leq \theta \leq 1$. This enables the application of Young's inequality above, within the required restriction (5.26), enforced by the GNS inequality when $n = 1, 2$.

Thus, we have

Let $\alpha = \frac{l}{2} < 1$. We have that $||v||^{2}_{2} \rightarrow 0$ as $t \rightarrow T^{*} < \infty$, for appropriately chosen initial data, in analogy with the ODE,

Since $L^{2}(\Omega)$ convergence implies uniform convergence on $\Omega$, which is closed and bounded, we see that for sufficiently chosen data $(u, v) \rightarrow (u^{*}(x), 0)$ uniformly, and this occurs in finite time. However, if $p = 1$, classical results via Theorem 5.7 ^[61], would imply the same data would have converged to $(0, v^{*}(x))$. This completes the proof. □

Remark 7. Theorem 5.8 tells us that in the spatially explicit setting with non-constant growth rate/resource distribution, once there is the slightest amount of flooding ($0 < q < 1$), there is an Allee threshold ($0 < p < 1$) associated with this flooding intensity, such that the population subject to the Allee effect will go extinct in finite time for certain initial populations.

Lemma 5.9. Consider (5.13) and (5.14), where $\Omega \subset \mathbb{R}^{n}$, $n = 1, 2$, is a bounded domain with smooth boundary, $0 < q < 1$. For positive initial data $(u_{0}(x), v_{0}(x))$, sufficiently small s.t $||v_{0}||^{2}_{2} < \left(q||m(x)||_{\infty}\right)^{-\frac{1}{1-p}}$, and a sufficiently chosen $p \in (0, 1)$, the solution $(u, v)$ to (5.13) and (5.14) will converge to $(u^{*}(x), 0)$ in finite time.

Proof. Consider the estimate for $v$,

Choosing, $y = ||v||^{2}_{2}$ yields

Where $C_{1} = 2||m(x)||_{\infty}$ and $C_{2} = 2$. We consider the ODE

Note via simple comparison, $y \leq V$, for all time $t$. Set $V = \frac{1}{U}$, then,

This yields,

Clearly, we see that $U$ will blow up in finite time for sufficiently large initial condition $U(0)$. That is if $U(0)$ is chosen s.t.,

Then

However, since $V = \frac{1}{U}$, we have,

That is, $V$ will go extinct in finite time, for initial data small enough, that is given by,

Since $y\leq V$, $v$ will go extinct in finite time as well. That is,

for initial condition chosen s.t.,

This implies setting $p = \alpha$ if

Then $v$ will go extinct in finite time, from which the convergence of $u$ to $u^{*}$ follows. This proves the lemma. □

6.   Application to the soybean aphid

6.1.   Background: The single biotype case

Kindlmann et al. ^[66] were among the first to propose a model for the population dynamics of Aphids via a system of differential equations,

$h(t)$ is the cumulative population density of a single aphid biotype at time $t$; $x(t)$ is the population density at time $t$. $a$ is a scalar constant, and $r$ is the growth/fitness rate of the aphids. The aphid population initially rises due to the linear growth term, but as the cumulative density becomes greater than the growth rate $r$, the population subsequently decreases due to the effects of competition. This results in a hump-shaped population density over time, typical of the boom-bust type scenarios witnessed with aphid populations ^[66].

Lemma 6.1. The equilibrium $(h^*, 0)$ for system (6.1) is non-hyperbolic.

Proof. The detailed proof can be found in Appendix (A.6) □

6.2.   Modeling virulent or avirulent aphids: Single biotypes on a resistant or susceptible plant

Soybean aphid biotypes are classified based on their ability to colonize soybean varieties expressing Rag, which is an acronym for "Resistance to Aphis glycines". For example, soybean aphid biotype 1 is susceptible to all Rag-genes; therefore, it is called avirulent. Biotype 2 is virulent to Rag1 ^[67], biotype 3 is virulent to Rag2 ^[68], and biotype 4 is virulent to both Rag1 and Rag2, capable of surviving on Rag1 and/or Rag2 plants, ^[69]. These four soybean aphid biotypes have been found throughout the soybean production areas in the Midwestern US ^[28,29]. Our aim in this section is to model the effect of excessive flooding, driven by climate change, on soybean aphid population dynamics. Herein, we aim to use the experimental results of Section 2.2. In these experiments, only a single bio-type (virulent or avirulent) is placed either on a susceptible or resistant plant - in both a control and a flooded scenario ^[34]. Recall what is found empirically is that there is a 28.4 $\%$ decrease in avirulent soybean aphid populations due to flooding on a susceptible soybean plant and a 44.1 $\%$ reduction on a resistant soybean plant. We aim to mimic this difference using our single species model (2.4), described earlier.

We first recap the classical result on the aphid population model 6.1,

Lemma 6.2. The equilibrium $(h^*, 0)$ for system (6.1) is globally attracting.

The proof of the above is given in ^[66].

We next consider (2.4), when $u = 0$, this yields, $\dfrac{dv}{dt} = v\left(a_{2}q - c(1-q) v^{p-1} - b_{2}v\right)$. Rewriting this yields,

We now model the above similar to (6.1), to enable boom-bust dynamics. For this the logistic term in (6.2), $\left(a_{2}q- b_{2}v\right) v$, is modified, to enable a boom-bust. To this end, we replace $b_{2}v$ with $b_{2}h$, where $h(t) = \int^{t}_{0}v(s)ds$, that is the cumulative density of $v$. In this setting the modified logistic term becomes $\left(a_{2}q- b_{2}h\right)v$. We also aim to reduce the number of parameters for tractability. The following model draws from the single biotype model ^[66] to model either a biotype 1 (avirulent) or a biotype 4 (virulent) aphid sub-population, attempting to colonize a susceptible or a resistant soybean host plant,

There are essentially only parameters $a, \alpha, p, r$, in the above. Here, $x = x_A(t)$, when we consider an avirulent aphid population density and $x = x_V(t)$, we consider a virulent aphid population density ^[32]. h is the cumulative population density of either avirulent or virulent aphids at time t; $r_1, r_2$ are the growth rates of the avirulent and virulent aphids, respectively, with typically $r_1 > r_2$ on a susceptible plant, and $r_1 < r_2$ on a resistant plant; a is a scaling constant relating prey cumulative density to its own dynamics. $0 \leq \alpha \leq 1$, is the flooding intensity. The case $\alpha = 0$ is the no flood (control) case, with $\alpha = 1$, the case where flood has maximum severity. We multiply the growth rate $r$ by $1-\alpha$ to model the adverse effect of flood on the growth rate of the aphid. $0 < p < 1$, models the strength of the FTE term. $p$ closer to 0 drives solutions down slower than $p$ closer to 1. This can be used to model the greater adverse effect on the avirulent populations in flooded scenarios than virulent ones, so $p_{1} > p_{2}$. Also note that if h is replaced by x in (6.3), then we recover (6.2), as the special case of (2.4), with $u = 0$, which is the basis of our modeling premise.

To this end, a theoretical result is stated,

Lemma 6.3. Consider the equilibrium $\left(h^*, 0\right)$ for system (6.3). There exists certain initial data that is attracted to this equilibrium in finite time.

Proof. The detailed proof can be found in Appendix (A.7) □

Next, based on recent empirical evidence ^[34], showing a clear fitness differential between virulent and avirulent biotypes due to flooding events, we run several numerical simulations.

6.3.   Numerical simulations

The parameters $r, p, \ \text{and} \ \alpha$ are chosen to replicate the dynamics observed in the experiments. For the first ten days, the population remains close to the numbers reported in the experiments (Figure 1). The model, however, runs for extended periods, so the variations in the peak population density in both control and flooding scenarios can be seen. Different values of $r$ are chosen for a specific biotype depending on whether the plant is susceptible or resistant. To simulate the flooding scenario, we apply $p$ and $\alpha$. However, in the control scenario, $\alpha$ is set to zero to eliminate any flooding effect. As mentioned in the results, the experiments were conducted for ten days, showing a $28.4 \%$ decrease in avirulent soybean aphid populations due to flooding on a susceptible soybean plant and a $44.1\%$ reduction on resistant soybean plant. From Figure 11(a), a $10.47\%$ decrease in the peak population of avirulent aphids is observed under flooding on a susceptible soybean plant when the model is simulated for 60 days. For the resistant soybean plant, the value of $r$ is chosen to mimic the experimental results; the model runs for 120 days as the peak occurs much later in the season. Figure 11(b) shows a decrease of $43.5 \%$ in the peak population density of avirulent aphids on a resistant soybean plant. Furthermore, Figure 12(a) and Figure 12(b) represent the dynamics of virulent aphids on susceptible and resistant soybean plants, respectively. The peak population decreases by $7.5 \%$ on a susceptible plant and by $11.01 \%$ on a resistant plant due to flooding.

7.   Discussion and conclusions

In this work, we propose a model for two competing species (or sub-populations) with differential fitness, driven in part by climatic changes - wherein one of the species is less fit than the other. This is modeled specifically via an Allee effect with finite time extinction, which affects the less fit sub-population. Several interesting results are derived. These are best summarized via Theorem 3.10, where based on parametric restrictions and the severity of the Allee term - that is, the exponent $p$ as well as the proportion of the population, $q$ that is affected differentially, one can get a complete change from classical Lotka-Volterra dynamics. One can avoid competitive exclusion and also have bi-stability dynamics; (Figure 9(b)). Also note that several bifurcation results are reported. Depending on the value of the exponent $p$ as well as the proportion of the population, when $q$ is differentially affected, one can have a saddle-node bifurcation in the interior (Figure 7(a)). Also, if one enables collision with boundary equilibrium points, a pitchfork bifurcation is possible; (Figure 8(b)). The finite time extinction via the Allee term is best quantified via Theorem 3.6, where an explicit relation for the smallness of the initial data is given to yield FTE. To this end, (Figure 5).

We also propose a proof of concept model, via (6.3), to mimic recent empirical evidence ^[34], where it is reported that avirulent biotypes may be (relatively) less fit on a flooded soybean host plant than their virulent counterparts. Furthermore, on a flooded susceptible host plant, avirulent aphid populations decrease up to 40 $\%$ from the no flood controls, whereas there is no significant decrease in virulent populations ^[34]. To model this scenario, our approach is to use the earlier defined parameters $\alpha, r, p$ in (6.3) to mimic these experiments. To reiterate, $r$ is the growth rate, $r_{1}$ is the growth rate of the avirulent, and $r_{2}$ is the virulent aphid, respectively, with typically $r_{1} > r_{2}$ on a susceptible plant, and $r_{1} < r_{2}$ on a resistant plant. $0 \leq \alpha \leq 1$ is the flooding intensity. $\alpha = 0$ is the no flood (control) case. With $\alpha = 1$, the maximum severity of the flood is possible. We multiply the growth rate $r$ by $1-\alpha$ to model the adverse effect of flood on the growth rate of the aphid as well. Also, $p$ models the degradation/differential fitness due to the effect of plant virulence on the different biotypes. This term causes population decay at the rate $\approx (T^{*} - t)^{\frac{1}{1-p}}$, where $T^{*}$ is the finite extinction time. Thus, $p > \frac{1}{2}$ (closer to 1) would model the scenario of avirulent aphids on resistant plant (including a flooded plant, as it is seen their populations drop severely due to flooding) and $p < \frac{1}{2}$ (closer to 0) would model the scenario of virulent aphids on a resistant plant (including a flooded plant, as it is seen their populations are less affected due to flooding). These simulations are consistent with empirical results from ^[34]. To reiterate, it is observed in experiments that avirulent soybean aphid populations decrease by 28.4 $\%$ with flooding on susceptible soybean and by 44.1$\%$ with flooding on Rag1+2 soybean. In contrast to these results, virulent aphids exhibited no significant response to flooding on either soybean variety (F$_{1, 32.2}$ = 0.275, P = 0.604; Figure 1). In the simulations, avirulent populations are seen to decrease by over 40 $\%$ on resistant plants (Figure 11(b)), whereas the virulent populations did not have any significant decrease (Figure 12(b)). Without loss of generality, we assume a few different growth rates $r$ in the biotypes in our numerical simulations. However, for the future, we will perform an exhaustive parameter sweep for different growth rates. All in all, the simulations via (6.3) are a stepping stone to guide future experimental work, where both aphid biotypes, virulent and avirulent, will be placed simultaneously on a single plant. Once these experiments are performed, further simulations via (2.5) will follow.

Moreover, other variable carrying capacity aphid population models can be explored, where the "variability" ties in closer to actual past flooding data and future projections. However, one must approach such models ^[70] with caution due to possible blow-up/explosive dynamics ^[71]. Note our analysis herein assumes a susceptible host soybean plant. Future work can include a resistant soybean plant, where virulent and avirulent aphids are both trying to colonize, and various dynamics are at play that cannot be described by earlier models. First, the virulent and avirulent are in direct competition for space, similar to inter-species competition. The virulent aphids are also in competition for space with other virulent aphids, as avirulent aphids are in competition for space with other avirulent aphids, similar to the intraspecies competition. These are direct effects of competition. Note that on a resistant plant, both the avirulent and virulent aphids are able to weaken the plant's defenses via feeding facilitation. However, for the avirulent aphid, this occurs only if it arrives in sufficiently large numbers ^[30]. Thus, there is a definite resistant level in the plant that is dependent on initial avirulent aphid density. If the avirulent aphids arrive in sufficient numbers above this level, they could colonize a resistant plant, but below this will die out–this is very similar to a resistance effect in ecology.

A more specific mechanism that inducts susceptibility on a resistant plant is the obviation of traits that confer resistance to the herbivore ^[30,72]. This mechanism requires a subset of the herbivore population that is virulent and capable of surviving on the resistant genotype of the host plant. By obviating the resistance through a physiological change to the plant, avirulent subpopulations can now survive on the resistant plant-this, then, is an indirect form of cooperation at play. Thus, the plants' resistance is a dynamic process ^[54], dependent on the presence and densities of these biotypes. Note that although such mechanisms are typical of a resistant host plant, this has not been considered in this work. However, modeling the differential fitness and subsequent variable carrying capacity effects on a resistant plant will make for interesting future work. To this end, various refuge concepts, such as the "within plant refuge", have been explored (albeit on a resistant host plant) ^[54]. In light of our results, one might consider how to build such refuges in the advent of climatic changes and higher precipitation, that is, flooding events. Clearly, flooding leads to a strong decrease in avirulent populations as compared to virulent populations. From a pest management perspective, this points to devising several strategies if one wants to keep the spread of virulence in check in the advent of climate change and increased flooding events. For example, an increase in the prevalence of virulent aphid biotypes under frequent flooding could highlight a need for increased scouting on aphid-resistant soybean varieties, with insecticide sprays necessary when populations increase economically damaging thresholds. Long-term, it may be possible to combat this increase in virulence by selectively breeding for simultaneous expression of aphid resistance and flood tolerance traits within soybean cultivars. Biological control from naturally occurring parasitoids and predators may also help to keep virulent aphid populations in check. Thus, considering the impact of natural enemies on aphid populations, and in particular, the dynamics of the plant-pest-natural enemy interactions in flooded scenarios vs not flooded ones, becomes important to study further. Dispersal is an important ecological feature that finds rich application in many areas ^[1,63,73]. Thus, the spatially explicit model for 2.5 is also considered for both the constant growth/resource case as well as the spatially inhomogeneous case. Here, we see that the effect of the FTE Allee term can lead to possible patterns (see Lemma 5.6). Also, in the spatially inhomogeneous case, one can obtain finite time extinction of the $v$ species, depending on the initial condition chosen (see Theorem 5.8 and Figure 10(d), (c)). These results tell us that depending on the choice of parameters, the spatially explicit model can predict results different than the ODE system (2.5)-in particular, pattern formation or "patchiness" of the populations is possible.

Current experimental work shows that control by parasitoids can keep pest populations lower than predator based control. Additionally, empirical observations show a numerical (but not statistically significant) increase in parasitoid attack rate under flooding conditions. Also, an increase in adult emergence rates with flooding (regardless of biotype) is also observed. Thus, it seems that flooding stress in soybeans is not going to affect biological control–via the use of parasitoids at least. Modeling these increased attack rates under flood conditions would make for interesting future work. In particular, it would be interesting to see if such a dynamic could balance the adverse effect of flooding on avirulent populations and what the overall dynamics of such mathematical models would predict. Researchers can also investigate models where a faction or group within the $v$ population decides to divert from the incentives of the group as a collective whole. Such competitive mechanisms have been under intense recent investigation ^[5]. Also, one can explore different formulations/motivations that would be more representative of such reaction terms. All in all, this work provides a segue into further empirical studies, where similar experiments, such those in ^[34], can be conducted, but with both biotypes on a single flooded host plant. The results therein would complete the feedback loop between empirical studies and modeling efforts so as to advise better pest management tactics for the soybean aphid in the event of ongoing global climatic changes.

Use of AI tools declaration

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

Acknowledgments

This work is supported by Agricultural and Food Research Initiative grant no. 2023-67013-39157 from the USDA National Institute of Food and Agriculture. AB is elated to acknowledge the motivating presence of Riti Sardar.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.1.   Proof of Lemma 3.1

Proof. Studying the boundary equilibrium $E_v(0, \bar{v})$, the nullcline (3.2) simplifies to $\phi(\bar{v})$ given by,

To study the roots of the polynomial $\phi(\bar{v})$, we study the monotonicity of the polynomial, which is described by,

Thus, the extrema is attainable when $\phi^{'}(\bar{v}) = 0$. We see, it is possible to attain extrema at $v = v_{\phi}$ where $v_{\phi}^{p-2} = \frac{b_2}{c(1-q)(1-p)}$ is positive as $b_2 > 0$ and $0 < p, q < 1$.

As $\phi^{''}(v_{\phi}) = c(1-q)(1-p)(p-2)v_\phi^{p-3} < 0$, we have a maxima at $v = v_\phi$. Now, as $v_\phi$ is unique for the fixed parameter set in the invariant set, we have unique maxima at $v = v_\phi$.

As there is one maxima, we can deduce there can be at most two roots of $\phi(\bar{v})$ and may have no roots depending on the functional value $\phi(v_{\phi})$. Computing we get, $\phi(v_{\phi}) = a_2q-b_2v_\phi(\frac{2-p}{1-p})$.

We know, $\lim_{v\to 0} \phi(v)$ = $\lim_{v\to \infty} \phi(v)$ = -$\infty$.

So, if $\phi(v_{\phi}) > 0$ i.e., $v_\phi < \frac{a_2q(1-p)}{b_2(2-p)}$ then there exists two boundary equilibrium points.

When $v_\phi > \frac{a_2q(1-p)}{b_2(2-p)}$, then there are no boundary equilibrium of the form $E_v(0, \bar{v})$. □

A.2.   Proof of Lemma 3.2

Proof. Let there exist a positive $v^*$ which is a root of the polynomial $\phi(v)$ given by,

To understand the monotonicity of the curve, we study the slope given by,

We know that, $0 < p, q < 1$ and $v^* > 0$ so $cb_1(p-1)(1-q)v^{*(p-2)} < 0$. Thus, $\phi'(v^*) > 0$ which shows $\phi(v^*)$ is a monotonically increasing curve. For $v\to 0$, $\phi(v) \to -\infty$ and for $v\to \infty$, $\phi(v) \to \infty$ so $\min(\phi(v)) < 0$. Thus, there exists a positive $z_1$ such that $\phi(z_1) < 0$ and $\exists$ positive $z_2$ such that $\phi(z_2) > 0$. By using the Intermediate Value theorem, we can conclude that there exists a $v^* \in (z_1, z_2)$ such that $\phi(v^*) = 0$ and as $\phi(v)$ is monotonically increasing so $v^*$ is unique. Thus, if $(b_1b_2-c_1c_2) \leq 0$, then there can only be one positive $v^*$ such that $\phi(v^*) = 0$. □

A.3.   Proof of Lemma 3.7

Proof. According to Lemma 3.1, we can have two, one, or no boundary equilibrium points of the form $E_v(0, \bar{v})$ depending on the functional value of the nullcline. We investigate when we have two boundary equilibrium points $E_{v_1}(0, \bar{v_1})$ and $E_{v_2}(0, \bar{v_2})$. Without loss of generality, let us assume $v_1 < v_2$. By Lemma 3.1, we know that $v_\phi$ is the maxima of the nullcline. By the Intermediate Value Theorem, we can assume $v_1 < v_\phi < v_2$.

To study the stability of the equilibrium points, we simplify the Jacobian matrix and have,

The eigenvalues for the system are $\lambda_1 = a_1-c_1v$ and $\lambda_2 = v(-b_2+c(1-p)(1-q)v^{(p-2)})$.

We study the case when there exist two boundary equilibrium points. When $v_2 > v_\phi$, where $v_{\phi}^{(p-2)} = \frac{b_2}{c(1-q)(1-p)}$, then,

$= v_2c(1-p)(1-q)(v_2^{p-2}-v_\phi^{p-2})$.

As $v_\phi < v_2$ and $p-2 < 0$ we have $\lambda_2 < 0$. Similarly, for $v_1 < v_{\phi}$ we have $\lambda_2 > 0$.

If we simplify $\lambda_1 = a_1-c_1v$, then $\lambda_1 < 0$ if $v > \frac{a_1}{c_1}$ for any equilibrium point $v$. □

A.4.   Proof of Lemma 3.8

Proof. According to Lemma 3.2 if $(b_1b_2-c_1c_2)\leq 0$ there exists an unique equilibrium $E^*(u^*, v^*)$. The simplified Jacobian matrix for the unique equilibrium is given by,

The determinant of the Jacobian matrix can be given by,

As $(b_1b_2-c_1c_2)\leq 0$ and $0 < p, q < 1$ so det$(J(E^*))$ is always negative.

Thus, $E(u^*, v^*)$ is always an unstable point. □

A.5.   Proof of Lemma 3.9

Proof. According to Theorem 3.4, if $(c_1c_2-b_1b_2) < 0$ and $(a_2b_1q-c_2a_1)+(c_1c_2-b_1b_2)(1+\frac{1}{1-p})v_{max} > 0$ where $v_{max} = -\frac{c(1-q)(1-p)}{(c_1c_2-b_1b_2)}^{\frac{1}{2-p}}$ then there exists two interior equilibrium points. We first study when we can get a stable interior equilibrium.

The conditions to be met for a stable equilibrium point are:

1) Existence condition : $(a_2b_1q-c_2a_1)+(c_1c_2-b_1b_2)(1+\frac{1}{1-p})v_{max} > 0$

2) Trace : $-b_1u^*-b_2v^*+c(1-p)(1-q)v^{*(p-1)} < 0$

3) Determinant : $(-b_1u^*)(-b_2v^*+c(1-p)(1-q)v^{*(p-1)})-c_1c_2u^*v^* > 0$

We simplify trace and determinant using the nullclines conditions as follows:

We start with simplifying the trace using (A.1),

Now, simplifying determinant,

Using (A.1) we get,

Combining the results (A.2) and (A.3), we can say that the interior equilibrium point is locally stable if,

From the existing condition, we get that, $q > \dfrac{1}{a_2b_1(1-p)}\{(b_1b_2-c_1c_2)(2-p)v_{max}+c_2a_1(1-p)\}$ □

A.6.   Proof of Lemma 6.1

Proof. From $\dfrac{dh}{dt} = 0 = f_1$ we have,

Using $\dfrac{dx}{dt} = 0 = f_2$ we have,

Therefore, we have a level set $(h^*, 0)$. Now, for the eigenvalues around this equilibrium point, we evaluate the corresponding Jacobian matrix,

The characteristic equation comes out as,

Both the eigenvalues are real, and one of them is always zero regardless of the sign of $(r-h^*)$. This proves the lemma. □

A.7.   Proof of Lemma 6.3

Proof. Consider the equation for $x$,

Set $x = \frac{1}{U}$, then,

We see that $U$ will blow up in finite time for sufficiently large initial condition $U(0)$. That is if $U(0)$ is chosen s.t., $-rU(0) +\alpha (U(0))^{(2-p)} \geq 0 \ \implies \alpha (U(0))^{(1-p)} \geq r$ then $\lim_{t \rightarrow T^{*} < \infty} U(t) = \infty$. However, since $x = \frac{1}{U}$, we have,

That is, $x$ will go extinct in finite time, for initial data small enough, that is given by,

This proves the lemma.□